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ABSTRACT 


In  this  thesis  the  problem  of  load  frequency  control  (LFC)  of  a 
single  area  and  an  interconnected  power  system  is  studied.  The  control 
strategy  which  has  been  chosen  here  adopts  what  is  implemented  in  the 
real  world  of  load -generation  control.  The  control  law  is  specified 
to  be  proportional  and/or  integral  of  the  area  control  error  (ACE) .  The 
area  control  error  is  a  measure  of  the  prevailing  generation  error. 
Different  strategies  have  been  proposed  to  handle  the  problem  of  LFC. 

The  LFC  problem  is  formulated  as  a  parameter  optimization  one,  in  order 
to  find  the  optimal  control  gains. 

Here,  the  tie-line  nonlinearity  has  been  taken  into  account.  In  this 
case,  the  system  nonlinearity  is  cast  into  the  quadratic  form  which 
facilitates  to  a  great  extent  the  convergence  of  the  numerical  solution 
of  the  parametric  optimization  problem. 

Also,  the  governor  dead  band  has  been  taken  into  consideration.  This 
problem  has  not  been  dealt  with  in  the  optimal  sense  until  this  work. 

The  computational  aspects  of  the  parametric  optimization  problem  and 
in  general  the  nonlinear  two-point  boundary  value  problem  is  discussed 
in  this  thesis. 
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Chapter  1 


The  Role  of  Load  Frequency  Control 
in  The  Power  System  Operation 


1.1  Introduction 

The  operation  of  an  electrical  power  system  may  be  viewed  as  a 
series  of  control  actions  taken  to  maintain  continuity  of  service  at 
standard  frequency  and  voltages  with  minimal  cost.  Automation  has  been 
always  integral  to  a  power  system  in  the  form  of  governor  action,  voltage 
regulation  and  protective  relaying  [1].  Further  automation  has  been 
achieved  in  the  form  of  automatic  generation  control,  load  frequency 
control  and  economic  dispatch. 

In  order  to  describe  the  load  frequency  control  role  in  the 
electric  power  system  operation,  one  may  use  the  concept  of  multilevel 
control  theory.  Accordingly,  the  system  mode  of  operation  may  be  decomposed 
into  the  following  [2]: 

(1)  Normal  mode 

(2)  Preventive  mode 

(3)  Emergency  mode 

(4)  Restorative  mode 

In  the  normal  mode  of  operation,  the  power  system  is  operated  so  that 
the  demands  of  all  the  customers  are  continuously  satisfied  at  standard 
frequency  and  voltages  with  minimal  cost.  Also,  the  system  spinning  reserve 
should  be  maintained  such  that  the  system  will  not  go  unstable  for  minor 
disturbances.  One  can  adopt  a  control  strategy  in  this  mode  to 
keep  the  voltage  at  approximately  standard  voltage 
keep  the  frequency  at  standard  frequency  ±  ef 
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keep  the  tie-line  power  interchange  at  its  scheduled  value 
meet  the  customer  demand  without  interruptions. 


2 


In  which  the  frequency  tolerance  can  be  taken  as 

ef  =  0.033  Hz 

according  to  the  North  American  Power  Systems  Interconnection  Committee 
(NAPSIC) . 

If  a  contingency  is  likely  to  occur  such  that  the  system  may  not  be 
able  to  return  to  the  normal  mode  of  operation,  then  the  system  is 
defined  to  be  in  the  preventive  mode  of  operation.  The  objective  of  the 
control  in  this  mode  is 

keep  the  frequency  at  60  ±  e^ 

-  keep  the  tie-line  power  flows  at  about  its  schedule  values 
keep  the  customer  demand  at  minimum  cost  subject  to  constraint 
on  the  amount  of  spinning  reserve. 

If  a  contingency  occurs  such  that  customer  demand  cannot  be 
maintained  at  prespecified  voltage  and  frequency,  since  otherwise  the 
electrical  system  will  lose  synchronism,  the  system  is  in  an  emergency 
mode  of  operation  and  the  control  objective  is 

-  to  keep  the  frequency  at  its  standard  value  ±  e^ 

to  maximize  supply  of  customer  demand  without  considering 
economic  factors. 

If  a  contingency  happens  and  service  to  some  customer  loads  is  lost, 
then  the  system  is  defined  to  be  in  the  restorative  operating  mode.  The 

control  objective  in  this  case  is 

to  bring  the  system  back  to  the  normal  or  preventive  mode  in 
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minimum  physically  possible  time. 

It  may  be  noted  that  in  all  modes  of  operation,  the  primary  goal 
of  the  control  strategy  is  to  try  to  keep  the  frequency  and  tie-line  power 
deviations  within  certain  tolerances. 

1.2  The  load  frequency  control  problem 

The  fact  that  electrical  energy  cannot  be  stored  in  large  quantities 
plays  an  important  role  as  far  as  power  generation  is  concerned.  Therefore, 
the  mechanical  output  of  system  turbines  must  be  continuously  controlled 
such  that  the  prevailing  electrical  generation  matches  exactly  the  consumer 
load  demand.  If  a  power  system  had  the  capability  to  do  so  instantaneously 
then  the  system  scheduled  frequency  and  scheduled  tie-line  power  flow  would 
remain  constant. 

Unfortunately,  such  a  capability  is  practically  impossible.  First, 
exact  prediction  of  load  demand  is  inconveivable .  Second;  it  is  not  feasible 
to  change  the  generated  power  to  match  load  demand  instantaneously  for 
mechanical  reasons,  e.g.,  turbine  output  rate  of  change  can  not  exceed  a 
certain  limit  due  to  thermal  stress  on  the  boilers.  Thus,  frequency  and 
scheduled  tie-line  power  deviations  occur  because  of  the  aforementioned 
reasons . 

1.3  The  need  for  load  frequency  control 

During  the  last  three  decades,  the  wide  spread  use  of  electric  clocks 
has  led  to  the  need  for  accurate  regulation  of  power  system  frequency.  There 
are  other  reasons  which  stem  from  the  operating  principle  of  power  systems 
in  addition.  A  major  reason  is  that  there  are  a  buy  a  sell  power  inter¬ 
change  agreements  between  adjacent  areas  to  supply  available  excess  lower 
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cost  energy  from  an  area  to  another  that  can  use  it  to  replace  its  own 
higher  cost  energy.  Exchange  of  power  between  adjacent  areas  of  a  power 
system  is  usually  governed  by  a  predetermined  schedule  so  that  during  a 
given  period  of  time  a  specified  amount  is  exchanged.  During  a  disturbance 
period,  the  outputs  of  all  system  generators  are  altered  since  each  area 
of  the  system  is  involved  in  accommodating  the  load  change.  Consequently, 
the  power  flowing  through  the  tie  lines  deviates  from  its  scheduled  values 
which  deviates  from  the  most  economic  operation  mode  of  the  system  [3]. 
Another  very  important  reason  for  keeping  the  frequency  within  a  certain 
tolerance  arises  from  the  operating  characteristic  of  the  system  turbines. 
When  a  system  operates  at  a  frequency  below  certain  limits,  some  types  of 
steam  turbines  undergo  excessive  vibrations  in  certain  turbine  rotor  stages 
which  result  in  a  metal  fatigue  and  blade  failure.  Moreover,  in  60Hz 
systems  sustained  operation  below  58.5Hz  of  such  turbines  may  limit  their 
operating  life  to  one  hour  at  full  load,  decreasing  progressively  with 
increasing  frequency  deviations  [A].  When  the  frequency  falls  below 
approximately  58.5Hz,  the  turbine  regulating  devices  will  be  fully  open 
and  the  generating  units  become  completely  loaded.  Further  decrease  in  the 
frequency  reduces  the  efficiency  of  the  auxiliary  mechanisms  at  steam  power 
stations,  especially  feed  pumps  whose  outputs  are  approximately  proportional 
to  one  third  of  the  frequency.  The  result  of  the  prolonged  operation  under 
this  condition  is  a  drop  in  the  output  generation  and  further  loss  of  power 
and  increase  in  frequency  deviation.  The  increase  in  frequency  deviation 
may  cause  an  avalanche  of  event  whose  nature  may  jeopardize  the  integrit\ 
of  operation  of  the  power  system. 
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1.4  Scope  of  investigation 

This  study  was  motivated  by  efforts  to  design  an  optimal  controller 
for  the  problem  of  load  frequency  control  [5-28]  of  electric  power  systems. 

In  chapter  3,  the  case  of  designing  an  optimal  controller  for  a 
single  area,  steam  and  hydro,  power  system  is  considered.  Different 
optimal  controllers  are  designed  in  a  systematic  way  so  as  to  reduce  the 
system  transients  and  bring  the  system  back  to  its  normal  state. 

In  this  case,  load  changes  are  accommodated  by  the  system  as  a  whole, 
regardless  of  where  on  the  system  they  may  occur.  The  control  strategy 
is  referred  to  as  Flat-Frequency-Control.  The  control  strategy  is  a 
proportional  and/or  integral  of  the  system  frequency  deviation.  It  is 
shown  that  the  optimal  control  gains  are  independent  of  the  load  variations. 

The  problem  of  LFC  involves  many  nonlinearities  in  the  system  dynamics. 

The  most  significant  of  these  arises  from  the  tie-lines  connecting  the  areas 
of  an  interconnected  system,  and  the  dead  band  of  the  area’s  governor.  The 
former  overides  when  the  system  is  subjected  to  large  disturbances  while 
the  latter  is  dominant  for  small  system  disturbances. 

In  chapter  4,  the  problem  of  linear  and  nonlinear  LFC  of  an  interconnected 
two  area  mixed  power  syster  is  studied. 

On  mode  modern  interconnected  power  systems,  the  basis  for  accomplishing 
the  regulating  responsibilities  is  the  Tie-Line-Bias  concept  (conventional 
LFC)  introduced  25  years  ago  [20].  When  operating  in  accordance  with  the 
principle  of  conventional  control,  each  area  of  an  interconnected  power 
system  attempts  to  regulate  its  ACE  to  zero.  This  can  be  achieved  by 
employing  a  proportional— plus-integral  control  of  the  area  control  error 
(ACE)  provided  that  each  area  can  fully  accommodate  its  load  changing.  The 
ACE  of  an  area  of  an  interconnected  power  system  is  a  measure  of  the 
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mismatch  of  the  area  output  generation  and  load  demand.  The  problem  of 
selecting  the  optimal  control  parameters  of  the  conventional  control 
strategy  was  studied  in  detail  by  Kirchmayer  [  14  ] .  In  [14],  the  problem 
of  a  two  area  interconnected  power  system  was  simulated  on  an  analog  computer. 
Then,  by  trial  and  error,  the  control  parameters  which  introduce  reasonable 
maximum  overshoot,  rise  time,  settling  time,  etc.,  for  the  frequency  and 
tie-line  power  deviations  were  considered  optimal. 

Recently,  several  attempts  have  been  made  to  employ  modern  control 
techniques  to  the  problem  of  load  frequency  control.  A  significant  group 
of  results  covering  various  approaches  to  design  an  optimal  controller  for 
the  problem  are  available  in  the  literature  [5-20].  Elgred  and  Fosha  [7] 
implemented  an  optimal  linear  regulator  theory  to  design  the  supplementary 
regulators  of  a  two  steam  area  interconnected  power  system. 

The  developed  control  was  a  function  of  all  the  systems  state  variables 
and  load  disturbances.  Therefore,  it  was  necessary  to  introduce  an 
observer  in  the  system  dynamics  to  realize  the  unmeasurable  states  and 
the  load  disturbances  [8].  Cavin,  et.  al  [44],  considered  the  LFC  problem 
from  the  stochastic  viewpoint.  In  [44],  the  system  was  assumed  to  be 
subjected  to  a  white  gaussian  random  variable  with  zero  mean  and  known 
covariance.  The  control  law  was  also  a  function  of  all  the  system’s  state 
variables,  however,  the  Kallman  filter  was  adopted  to  estimate  the  unobservable 
states.  The  assumption  that  the  distribution  of  load  disturbance  is  white 
noise  is  far  from  being  adequate  in  describing  disturbances  found  in 
electric  power  systems.  Also,  the  determination  of  a  meaningful  contriance 
matrix  in  practice  is  very  difficult  as  detailed  statistical  data  about  a 
plant  and  its  measured  noise  is  required  and  such  data  is  generally  not 
available  for  a  power  system.  Although  a  power  system  load  disturbances 
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are,  in  general,  random  by  nature,  it  is  reasonable  to  assume  that  a  power 
system  load  disturbances  are  deterministic  as  modern  LFC  systems  are  designed 
with  filters  [5]  in  order  to  remove  the  purely  random  portion  of  the  area's 
load  disturbances  leaving  the  deterministic  components  for  the  system  units 
to  regulate. 

So  for  all  the  work  mentioned  before  have  attempted  to  control  the 
frequency  and  the  tie-line  power  deviations  independently  rather  than 
regulating  the  area  control  error  (ACE)  based  on  the  frequency  bias  which 
is  the  goal  of  tie-line-bias  control  used  in  practice.  Calovic  et.  al . , 
[10-11]  and  [51]  tried  to  narrow  the  gap  between  the  conventional  and  advanced 
LFC  by  letting  the  integral  of  the  ACE  to  be  a  part  of  the  control  strategy 
[12],  In  [10],  the  proportional  part  was  a  function  of  all  the  system 
state  variables.  In  [51],  both  the  proportional  and  the  intergap  part  was 
a  function  of  the  output  state  variables.  The  necessary  conditions  to 
validate  the  developed  control  law  depends  upon  the  choice  of  a  certain 
weighting  matrix.  However,  as  reported  by  the  authors,  there  is  no 
systematic  way  of  computing  this  matrix.  The  economic  control  is  a  "tertiary" 
control,  i.e.,  relatively  slow  and  is  based  on  quasi-static  description  of 
generation  and  demand  [70].  Therefore,  the  economic  dispatch  should  be 
treated  separately  from  load  frequency  control. 

In  this  chapter,  the  tie-bias-control  concept  which  is  extensively  used 
in  practice  is  adopted  as  a  part  of  the  proposed  control  strategy.  The 
control  law  of  each  area  of  a  power  system  is  assumed  to  be  proportional- 
plus-integral  of  the  area  control  error  (ACE).  Moreover,  the  strategy 
adopts  the  principle  that  only  the  disturbed  area  supplementary  regulator 
will  react  in  response  to  its  area  load  change  (nonintervention  principle); 
provided  that  the  disturbed  area  has  the  capability  to  accommodate  its  load 
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change,  is  adopted  here  for  the  problem  of  LFC  of  interconnected  power 
system. 

Here  the  problem  is  formulated  as  a  parameter  optimization  one.  It 
has  been  found  that  the  optimal  control  parameters  are  independent  of  the 
load  variations  in  contrast  to  [7],  and  therefore  there  is  no  need  to 
identify  it  [8].  Moreover,  the  proposed  strategy  adopts  the  principle  that 
only  the  disturbed  area  supplementary  regulator  will  react  in  response  to 
its  area  load  change  (non-intervention  principle) ;  provided  that  the 
disturbed  area  has  the  capability  to  accommodate  its  load  change.  This  is 
in  contrast  to  most  of  the  recent  proposed  strategies  [5-28]. 

Here,  the  problem  of  LFC  of  interconnected  power  system;  either  linear 
or  nonlinear,  is  formulated  as  a  parameter  optimization  one.  In  the  case 
of  the  linear  one,  the  optimal  control  parameters  have  been  found  to  be 
independent  of  the  load  variations  and  therefore  there  is  no  need  for 
incorporating  an  identifier  (observer)  in  the  system  dynamics  to  estimate 
the  load  demand  [8]. 

Many  efforts  have  been  done  to  study  the  nonlinear  version  of  the 
LFC  problem;  tie-line  nonlinearity.  Miniesy  and  Bohn  [16]  have  proposed  a 
two  level  controller  for  the  problem  of  nonlinear  load  frequency  control. 

The  first  level  is  a  local  feedback  control  for  each  area  of  the  interconnected 
system.  This  control  is  a  function  of  its  own  area  state  variables.  The 
second  level  is  an  intervention  open  loop  which  is  utilized  to  compensate 
for  neglecting  the  coupling  state  variable  between  the  areas  of  the  inter¬ 
connecting  system,  and  the  nonlinearities  due  to  the  tie-lines.  The  idea 
of  employing  a  multi-level  control  strategy  in  the  area  of  load  frequency 
control  is  not  adequate.  First,  it  increases  the  controller  design  complexity. 
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Second,  it  is  not  used  in  practice.  Doraiswami  [17]  has  considered  the 
problem  of  LFC  with  a  tie-line  nonlinearity  from  the  stochastic  viewpoint. 

An  observer  was  employed  to  implement  the  control  law. 

Here,  the  system  nonlinearity  is  cast  into  the  quadratic  form  and  the 
problem  is  formulated  as  a  parameter  optimization  one.  The  optimal  control 
parameters  have  been  found  to  be  dependent  on  the  load  variations.  However, 
it  is  shown  that  the  suboptimal  control  parameters  obtained  from  the  linear 
model,  are  good  enough  to  design  supplementary  area  regulators  with 
insignificant  (in  a  practical  sense)  degradation  of  system  performance. 

This  eliminates  the  necessity  of  designing  an  adaptive  controller  [17]  and 
therefore  simplifies  the  design. 

In  chapter  5,  the  governor  deadbank  of  a  single  area  steam  power 
system  is  studied. 

With  the  recognition  that  a  well-developed  linear  control  theory  now 
exists,  more  research  is  being  directed  towards  nonlinear  aspects  of  general 
control  systems.  Linear  control  theory  (classical  or  optimal)  suffers  from 
the  fundamental  criticism  that,  in  reality,  dynamical  systems  are  frequently 
subject  to  several  complicating  factors  which  may  invalidate  or  at  least 
severely  limit  its  applicability.  There  are  no  general  methods  for 
analysis  and  synthesis  of  nonlinear  control  systems  in  classical  or  modern 
control  theory.  For  example,  dealing  with  systems  with  control  signal 
saturation  (relay  control  systems  [71]  and  [72]  employing  modern  control 
theory  is  quite  different  from  dealing  with  systems  with,  e.g.,  output  feed¬ 
back  which  contains  a  functional  nonlinear  element  [73]. 

Here  we  deal  with  a  system  whose  dynamics  comprise  a  hysteresis  element. 
The  approach  developed  here  (chapter  5)  is  general  and  is  applied  to  the 
problem  of  finding  the  optimal  control  parameters  of  the  supplementary 
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regulator  of  single  area  steam  power  systems. 

The  problem  of  finding  the  optimal  control  parameters  has  been  reduced 
to  the  problem  of  solving  a  nonlinear  two-point-boundary-value  problem 
(TPBVP)  with  multiple  corners.  This  has  been  handled  by  employing  the 
Weierstrass-Erdmann  conditions  for  optimality  [29].  A  technique  is  developed 
to  detect  when  a  corner  and  hence  the  problem  can  be  treated  as  ordinary 
TPBVP . 


Chapter  2 


Parametric  Optimization 


2.1  Introduction 

Parametric  optimization  techniques  represent  a  wide  class  of  aids 
in  designing  control  systems.  The  tools  which  are  required  to  perform 
parameter  optimization  processes  depend,  among  other  things,  upon  the 
nature  of  the  input  signal  to  the  dynamic  system  under  consideration.  The 
nature  of  the  input  signal  to  any  dynamic  system  is  either  stochastic  or 
deterministic  processes.  In  this  dissertation,  systems  subjected  to 
deterministic  inputs  are  considered. 

2.2  The  problem  of  deterministic  parameter  optimization 

The  parameter  optimization  problem  considered  here  may  be  formulated 
in  terms  of  ordinary  differential  equations  of  the  form 


X(t)  =  f(X(t),C),  te[to,  tf] 


(2.1) 


Here,  t  is  the  independent  variable,  X  is  an  n-dimensional  vector  whose 
components  are  the  dependent  variables,  anc  C  is  an  m-dimensional  vector 
whose  components  represent  the  salient  variables  of  the  system.  C-vector 
is  assumed  to  be  time-independent,  that  is, 


C  =  0, 


(2.2) 


When  the  parameters  in  C  and  the  complete  initial  conditions 


(2.3) 
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are  known,  a  numerical  integration  of  Eqn.  (2.1)  produces  the  solution  X(t) 
on  the  time  interval  t^  <_  t  <_  t^.  This  initial  value  problem  can  readily 
be  solved  with  a  digital  or  analog  computer. 

On  the  other  hand,  the  parameter  optimization  problem  is  the  problem 
of  finding  the  parameter  vector  C  which  maximizes  or  minimizes  some  scalar 
performance  index 

t 

J  =  /  f  4>[X(t),  C , t ] dt  +  G[X(tf),C]  (2. A) 

t0 

subject  to  system  dynamic  constraints  given  by  Eqns .  (2.1)  -  (2.2).  It 
is  assumed  that  the  integrand  $  has  continuous  first  partial  derivatives 
with  respect  to  all  of  its  arguments;  t^  and  t^  are  fixed.  It  is  also 
assumed  that  the  trajectory  of  X(t)  have  only  piecewise-continuous  first 
derivative;  that  is,  X(t)  will  be  continuous  except  at  a  finite  number  of 
times  in  the  interval  (t^,t^). 

It  can  be  seen  that  Eqn.  (2.4)  is  sufficiently  general  to  investigate 
a  wide  class  of  practical  problems. 

Only  in  some  simple  parameteric  optimization  cases  is  it  possible  to 
obtain  an  analytic  solution,  but  in  most  cases  a  numerical  method  must  be 
used.  There  are  many  numerical  methods  available  for  this  purpose,  e.g. 
[30-32];  some  of  these  references  require  the  evaluation  of  the  performance 
index  only,  while  others  require  the  evaluation  of  the  gradient,  and 
possibly,  of  the  higher  order  derivatives  of  the  performance  index. 

The  method  of  gradients  which  has  been  applied  with  considerable 
successes  to  a  variety  of  problems,  is  adopted  here  to  solve  the  f orementioned 
parametric  optimization  problem.  The  method  of  gradient  requires 
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1  -  The  evaluation  of  the  performance  index,  and 

2  -  The  evaluation  of  the  gradient  of  the  performance 

index  with  respect  to  the  parameters. 

The  first  requirement  can  be  carried  out  by  using  a  numerical 
integration  process  by  augmenting  the  system  dynamic  equations,  Eq.  (2.1), 
by  adding  the  equation 


xq  =  *[X(t),C,t],  x0(tQ)  =  0 


(2.5) 


Then,  Eqns.  (2.1)  and  (2.4)  are  evaluated  in  the  forward-time  to  yield 
the  required  value  of  the  performance  index 


J  =  xQ(tf)  +  G[X(tf)C] 


(2.6) 


The  second  requirement,  the  evaluation  of  the  gradient  of  the 
performnce  index, 


J 

C  9C 


(2.7) 


can  be  carried  out  conveniently  in  terms  of  a  Hamiltonian  defined  by 


H  ■  t  +  A  f 


(2.8) 


In  Eqn .  (2.8),  A  is  referred  to  as  the  co-state  vector,  and,  is  defined  by 


.  SH  ,  .  N  3G  |  _  r 

A  ax  *  A(tf)  ax  'tf  Gxf 


3H 

ax  ’ 


(2.9) 


Then  the  gradient  of  the  performance  index  can  be  shown  to  be  given  by 
the  following 


J 


C 


t 


(2.10) 


0 


where 


(2.11) 


=  I 

C  9C  1  t 


f 


(2.12) 


and  it  can  be  evaluated  conveniently  by  introducing  the  new  variables 


Q  -  -  Hc-  Q<V  ■  Gc 


(2.13) 


Then,  the  integrating  of  Eqns .  (2.9)  and  (2.13)  backward  in  time,  yield 


(2.14) 


Jc  ' 


provided  that  the  state  vector  from  the  forward-time  integration  has  been 
stored.  The  scheme  of  forward-time  and  backward-time  integration  used  her 
is  typical  of  the  two-point  boundary  value  problems,  and  generally  arises 
in  the  optimization  problems  defined  on  finite  time  intervals. 

The  two-point  boundary  value  problem  arising  in  the  parametric 
optimization  problem  considered  here  can  be  described  as  follows: 


(i)  Forward-time  equations: 
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X  =  f[X(t),  C],  X(t0)  =  x0 

xQ=  4>[X(t),  C,  t],  x(tQ)  =  0 

The  value  of  the  performance  index 


(2.15) 


J  =  xQ(tf)  +  G[X(tf) ,  C] 


(2.16) 


can  be  found  by  solving  equations  (2.15)  forward  in  time. 


(ii)  Backward-time  equations: 


A  =  -  H„,  A(t,)  =  G 

A  f  X 

Q  =  "  V  Q<tf)  -  Gc 


(2.17) 


The  value  of  the  gradient  of  the  performance  index  with  respect  to  C 


Jc  «  Q(t0)  (2.18) 

can  be  found  by  solving  equations  (2.17)  backward  in  time. 

Once  the  gradient  of  the  performance  index  with  respect  to  the  vector 
C  is  found,  then  one  can  solve  the  two-point  boundary  value  problem 
described  by  Eqns .  (2.15)  and  (2.17),  by  using  successively  the  expression 


C  =  a  D  +  E 


(2.19) 
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where  D  is  a  feasible  direction  depends  on  J  ,  a  is  an  optimal  or 
subopt imal  step  size,  and  E  is  the  previous  estimate  of  the  vector  C. 

One  of  the  well  known  methods  to  determine  a  feasible  direction  D  is 
the  "Method  of  steepest  descent"  or  the  "method  of  gradient".  Recall  that 
the  direction  of  steepest  ascent  of  J  at  a  given  value  E  of  C  is  equal  to 
^E *  The  vector  “  therefore  points  in  the  direction  of  steepest  descent. 
Starting  with  an  estimate  parameter  vector  E  we  minimize  J  successively 
along  lines  in  directions  of  the  steepest  descent.  The  method  of  gradient 
is  useful  for  a  large  class  of  well  conditioned  problems.  However, 
experience  has  shown  that  the  method  can  be  extremely  slow  [ 33  ] .  Fortunately, 
there  is  a  nice  modification  of  the  gradient  method,  called  the  conjugate- 
gradient  method,  which  yields  the  optimum  solution  much  faster  than  the 
method  of  steepest  descent.  In  our  work  here,  the  "method  of  conjugate- 
gradient  is  employed  to  find  the  direction  D  appearing  in  Eqn.  (2.19). 

Further,  we  determine  the  optimal  step  size  a  by  solving  the  one¬ 
dimensional  minimization  problem 

min  J(E  +  aD)  (2.20) 

a  >  0 

2.3  Conjugate-gradient  algorithm 

The  method  of  conjugate-gradient  is  a  special  case  of  the  method  of 
conjugate  directions  where  the  successive  feasible  directions  are  related 
to  each  other  by  the  relationships 

dT  Q  D  .  =  0  ,  i, j  =  1,2,...,  i^j 


(2.21) 
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where  Q  is  any  positive  definite  matrix.  In  this  case  we  say  that  D 

i 

and  are  Q-conjugate.  The  method  of  conjugate-gradient  associates 
conjuancy  properties,  Eqn.  (2.21),  with  the  method  of  steepest  descent  in 
one  algorithm  in  order  to  achieve  both  efficency  and  reliability.  In  the 
case  of  optimizing  a  quadractic  performance  index,  the  method  of  conjugate- 
gradient  is  an  exact  one-dimensional  search  method  [33];  i.e.  it  does  not 
require  a  one— dimensional  search  in  conjunction  with  the  gradient— conjugate 
algorithm.  Note  that,  this  method  is  not  an  exact  one-dimensional  search 
method  for  a  nonquadratic  performance  index. 

In  the  case  of  nonquadratic  performance  index,  there  are  a  number  of 
algorithms  that  might  be  adopted  in  conjunction  with  the  conjugate-gradient 
algorithm.  Since  the  gradient  of  the  performance  index  has  to  be  found 
to  determine  the  feasible  directions,  it  is  reasonable  to  use  a  version 
of  a  one-dimensional  search  algorithm  based  on  the  derivatives  of  the 
performance  index. 

Here,  the  restrictive  extrapolation  and  interpolation  technique 
developed  by  Fletcher-Reeves  [33]  is  employed  to  find  the  optimal  step 
size  a. 

Let  vector  G  be  defined  as  the  performance  index  gradient  with 
respect  to  C;  that  is, 


i 


(2.22) 


Then  the  feasibility  direction  vectors  can  be  expressed  as 


D.  =  G.  +  {G.  G./G 


T  T 


i  i  i  i  i 


(2.23) 
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We  shall  refer  each  computation  for  a  feasible  direction  and  corresponding 
optimization  processes  in  one— dimension  as  an  iteration. 

Further  each  minimization  of  J(C  +  aD)  with  respect  to  a  >  0  will  be 
referred  to  as  a  trial.  The  general  trial  procedure  can  be  specialized 
by  selecting  initialization,  extrapolation,  and  interpolation  procedures. 

The  basic  steps  of  finding  the  optimal  parameters  by  the  method 
of  conjugate-gradient  are: 

(1)  Select  a  reasonable  value  of  the  parameter  C1,  and  store  it  in 
the  memory  of  the  digital  computer.  Let  the  iteration  index  i  be  zero. 

(2)  Using  the  nominal  parameter  C1,  integrate  the  state  equations 

(2.8)  from  tQ  to  tf  with  initial  conditions  X(tQ)  =  XQ  and  store  the 
resulting  trajectory  as  a  piecewise-constant  vector  function. 

(3)  Calculate  A  (t^)  and  Q1(t^).  Then,  integrate  equations  (2.15) 
backward  in  time  from  t^  to  t^  using  the  piecewise-constant  vector  stored 
in  step  (2).  Then,  evaluate  the  performance  index  gradient  with  respect 
to  parameter  vector  from  equation  (2.16)  and  store  the  vector  J  . 

(A)  If 


< 


£ 


1 


where 


is  a  prespecified  positive  constant, 


and 


2 


/  f  JT.  dt 

t  c1 

to  c 


(2.23) 


terminate  the  iterative  procedure.  If  the  iterative  stopping  criterion, 
equation  (2.22)  is  not  satisfied,  perform  the  trial  process: 

a  -  Select  a  reasonable  initial  value  of  the  step  size  aQ.  Here, 


* 
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quadratic  extrapolation  based  on  an  estimate  of  the  minimum  value  of  the 

performance  index;  J  .  ,  or  a  prespecified  a  .  ,  results  in  the  initialization 

mm  r  mm 


procedure  [32].  Take  aQ  =  0,  then  calculate 


J  .  -  JCC1) 

.  r  „  mm  , 

=  mm {2  - - -  ,  a__._  } 


T 

D  G. 

l 


mm 


(2.24) 


b  -  Using  the  step  size  as  in  equation  (2.24), 
This  is  done  through  steps  (2)  and  (3) . 
c  -  If 


calculate  J  . 

C^+ad 


1  e2  (2.25) 

where  is  a  given  positive  constant  <  1  terminate  the  trial  process 
and  return  to  step  (2).  If  the  trial  stopping  criterion,  equation  (2.25), 
is  not  satisfied,  a  restrictive  extrapolation  or  interpolation  is  implemented 
to  generate  a  new  step  size. 


T 

D  G  . 

cVt.D. 

_ 1  l 

T 

D  G  . 

c1 


a  =  Qq  +  (a^  -  a^)r  (2.26) 


In  the  case  of  J  (a, )<J  (a  )  <0,  the  performance  index  J(a)  is  approximated 

a  0  a  1  — 

by  a  quadratic  function  H(a),  of  the  form  [71] 


2 

H(a)  =  a  +  ba  +  ca  (2.27) 


If  the  values  of  J(aQ),  J  (aQ),  and  Ja(a1)  are  available,  then  the  parameters 
a,  b,  and  c  of  H  can  be  determined  by  solving  the  equations 


' 


i  t 
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2 

a  +  b  ao  +  c  ao  =  J^a0') 

b  +  2  c  aQ  =  Ja(a0)  (2.28) 

b  +  2  c  a.  =  J  (a. ) 

1  a  1 

The  extrapolating  quadratic  polynomial  can  be  written  accordingly  by 


H(a)  =  J(a  )  +  J  (a  )  [a-a  ]  + 
U  a  0  0 


0.5  (a-a0) 


[J  (a  )  -  J  (a  )] 
l  a  0 _ a  1 

ao  "  ai 


Condition  that  H  (a)  =  0,  yields 

a 


a  =  aQ  +  (a^  -  aQ)r 


(2.29) 


(2.30) 


where 


r  = 


a  0 


J  (a  )  -  J  (a  ) 
a  U  a  1 


(2.31) 


However,  if  J  (a„)  >  J  (a,)  then  a  does  not  exist,  and  if  J  (a  )  -  J  (a ,  ) 
a  0  —  a  1  aO  al 

is  small,  then  a  may  be  unacceptably  large.  For  these  reasons  we  restrict 
a  so  that  in  Eqn.  (2.26)  the  quantity  r  will  be  in  the  form 


r= 


f 


J  (O  i 

- - - — - - r-  for  J  (a  )  >  (1-  - )  J  (a  ) 

W  '  J„(Bl)  “  1  "  "max  °  0 


(2.32) 


“  '  — )  W 


max 


for  J  (a. ) 
a  1 


< 


max 
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and 


re (1 ,  r  ), 
max 


(2.33) 


where  r  is  chosen  as  follows 
max 


r  e(5,  10). 
max 


(2.34) 


In  the  case  of  the  conditions  that  J  (a_)<  0  and  J  (a, )>  0  are  satisfied 

a  0  a  1 

a  cubic  interpolation  can  be  implemented  to  find  the  optimal  step  size.  As 
before  J  (c  +  ad)  =  J(a)  can  be  approximated  by  a  cubic  function  H(a) 

[71]  in  the  form, 


2  3 

H(a)  =  a  +  ba  +  ca  +  da 


(2.35) 


Assuming  that  the  values  of  J(a  ),  J  (a  ) ,  J  (a  )  ,  and  J  (a  )  are  available, 

(J  cx  vj  -L  a  i 

then  the  parameters  of  H  can  be  determined  from  the  equations 


a  +  bJ (Uq)  +  c[J(aQ)]2  +  d[J(aQ)]3  =  J (aQ) 


a  +  bJ(a^  +  c[J(a^)]2  +  d[J(a^)]  -  J(a^) 


(2.36) 


b  +  2c  J (aQ)  +  3d[J(aQ)]  =  Ja(a0.) 


b  +  2c  J(a1)  +  3d[J(a1)]  =  Ja(a1> 


In  this  case,  the  interpolating  polynomial  H(a)  is  of  the  form 
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H(a)  =  J(a0)  +  VV^rV  '  [Z  +  Jo(a0)1(a‘a0) 


+  3  [Ja(“0)  +  J0(a)  +  2zHc<-a0) 


(2.37) 


H  (a)  =  0  yields 

a 


a  =  aQ  +  (a  -  0Lq)t 


(2.38) 


where 


r  = 


W  +  z  +  q 

W  +  W  +  2Z 


(2.39) 


in  which 


J  (a  )  -  J  (a.  ) 

Z  =  3  ^ -  +  J  (a  )  +  J  (a  ) 

al  ~  a0  a  0  a  1 


(2.40) 


2  1/2 
q  =  [Z  -  J  (a  )  J  (a  )]±/Z 
a  0  a  1 


(2.41) 


when 


J  (O  +  J  (a.  )  +  2Z  =  0 
a  0  a  1 


(2.42) 


then,  r  is  given  by 


r  = 


J  (°n> 
a  0 


W  -  J«(al) 


(2.43) 


23 


Therefore,  condition  H^(a)  —  0  yields  an  estimate  a  in  the  form  of 
(2.20)  with  r  as  follows 


r  =' 


W  +  z  +  q 

W  +  W  +  2z 


Ja(a0)  •  W 


for  Ja(aQ)  +  Ja(a1)  +  2z^  0 

for  J  (a  +  J  (a  )  +  2z  =  0 
a  0  a  1 


(2.44) 


d  -  Replace  by  and  ot^  by  a  and  return  to  step  b  in  the  trial 
process . 

2.4  Integration  methods 

It  has  been  shown  here,  that  each  computation  of  J(C)  and  requires  the 

computation  of  the  state  and  co-state  equations.  For  an  efficient 
computational  procedure,  then,  it  is  desirable  to  minimize  the  time  required 
for  the  integration  process.  This  reduction  can  be  achieved  by  employing 
a  predictor-corrector  type  of  integration  scheme. 

Predictor-corrector  methods  require  the  computation  of  the  derivatives 
two  times  in  each  integration  step,  as  compared  to  four  times  for  the 
standard  Rune-Kutta  method  [34].  Here  Hamming  predictor-corrector  scheme 
has  been  employed  to  perform  the  required  integration  processes  in  section 
(2.3). 

2 . 5  Remarks 

It  is  obvious  that  the  success  of  the  gradient  methods  of  optimization  is  somewhat 
dependent  on  judgement  and  intuition  of  the  user.  It  has  been  successfully 
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applied  to  some  practical  problems  such  as  the  problem  of  load  frequency 
control  which  is  the  subject  of  this  dissertation. 

A  computer  program  has  been  written  in  Fortran  IV  to  simulate  the 
proceeding  algorithm.  A  complete  list  of  the  code  is  given  in  Appendix  A. 
This  code  can  be  used  for  large  classes  of  optimal  control  design  problems. 


Chapter  3 

Load  Frequency  Control  of  Single  Area  Systems 

3.1  Introduction 

A  single  area  system  is  one  in  which  load  changes  are  accommodated  by 
the  system  as  a  whole,  regardless  of  where  they  occur  on  the  system.  Load 
changes  that  occur  in  any  part  of  the  system  may  be  absorbed  elsewhere 
within  the  system,  in  accordance  with  allocation  practices  prevailing 
at  that  particular  time.  No  one  part  of  the  system  is  expected  to  adjust 
its  own  generation  to  counteract  its  own  load  changes. 

Most  large  power  systems  are  geographically  wide-spread,  made  up  of 
different  operating  companies.  Experience  in  automatic  operation  of  these 
systems  sometimes  shows  a  need  to  assign  regulation  of  system 
frequency  to  one  of  these  companies  in  the  group  centrally  located  with 
respect  to  other  companies.  Regulation  of  the  remaining  companies  in  the 
system  depends  on  the  frequency  and  the  power  flow  of  the  tie-lines  connecting 
the  respective  companies  to  the  central  area.  Therefore,  the  central  area 
tie-line  power  flows  are  neither  scheduled  nor  controlled.  A  load  frequency 
control  which  adopts  complete  relaxation  of  the  tie-line  power  interchange 
between  the  companies  of  a  system  is  referred  to  as  Flat-Frequency  Control. 

This  chapter  is  devoted  to  the  study  of  load  frequency 
control  for  systems  or  areas  which  employ  a  flat  frequency  control  strategy. 
The  parameter  optimization  technique  developed  in  chapter  2 
is  employed  to  design  an  optimal  controller  for  the  frequency  control 
of  electric  power  generation. 
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3.2  Load  frequency  control  of  single  area  steam  power  system 
3.2.1  Problem  formulation 

A  typical  linearized  model  of  a  single  area  steam  power  system; 
SASPS,  is  shown  in  fig.  (1)  [15].  On  state  variable  form,  the  system 
dynamics  can  be  written  as  follows: 


X.  = 


A1X1  +  B  u  +  r 


Y!  “  c!x! 


where 


X  =  (Aw,  Ax  ,  Ap  ) 

g  g 


in  which 


(3.1) 


Aw  =  angular  frequency  deviation  rad/sec. 


Ax  =  deviation  in  governor  position  in  pu  power 
g 

Ap  =  deviation  in  turbine  output  power  in  pu  power 
g 

u  =  speed  changer  position  in  pu  power 


AL  =  step  load  change 


and 


-  G 

M 

0 

1 

M 

E 

1_ 

0 

T 

T 

g 

g 

1 

1 

T 

T\ 

t 

t 

l-Ue- 


AL 

1 

1 

K 

1 

TP  +  1 
g 

TP  +  1 

J - 

MP  +  G 

Governor 

Turbine 

Alternator 

Fig.  (1)  Block  diagram  of  a  single  steam  area  power  system 
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^=(100),  B1  =  (1,  0,  0),  and  T  =  (1,  0,  0) 


The  control  is  assumed  to  be  integral  or  proportional  -plus- 
integral  of  the  angular  frequency  deviation.  In  terms  of  the  system  state 
variables ,  the  control  law  can  be  written  as  follows t 

tf 

u  ■  -  kpxrki  /  x!dt  o.2) 

t0 


in  which 


=  proportional  gain 
kj  =  integral  gain 

the  control  parameters  k  and  k_  can  be  expressed  mathematically  as 

pi 

follows : 


K  =  0 


(3.3) 


where 


The  problem  posed  is  to  find  the  optimal  control  u  or  equivalently 

the  optimal  integral  gain  k^  and/or  the  proportional-plus-integral  control 

parameters  k  and  k,  which  minimize  the  cost  functional 
P  I 


1 


(xlQlX!  +  y  u2)dt 


(3.4) 


subject  to  the  dynamic  constraints  given  by  (3.1)  as  well  as  to  satisfy 


. 
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the  system  transient,  steady  state  and  stability  specifications.  Here, 

Q1  Js  a  posltlve  definite  diagonal  matrix,  y  is  equal  to  1.  By  adjusting  the 
element  values  of  we  can  weigh  the  relative  importance  of  each  of  the 
system  states.  Thus,  by  increasing  q^  we  attach  more  significance  to  x . , 
by  making  zero  we  indicate  the  state  x^  is  of  no  concern  whatsoever. 

Here  all  the  states  have  the  same  importance.  Therefore  Q  has  been  choosen 
to  be  the  identity  matrix. 

The  transient  requirement  is  met  when  the  generation  rate  limit  is 
kept  within  ¥  pu  power/sec.  The  generation  rate  limit  is  not  fixed  and 
varies  from  unit  to  unit.  The  transient  specification  can  be  mathematically 
expressed  in  terms  of  the  state  variables  as  follows: 


1 


(3.5) 


In  fact  the  R.H.S.  of  (3.5)  is  the  R.H.S.  of  the  third  equation  in  (3.1). 

The  steady  state  requirement  is  met  when  the  frequency  deviation  is 
equal  to  zero  at  the  end  of  the  control  period  provided  that  the  disturbed  area 
can  fully  accommodate  its  load  changes.  This  is  achieved  as  a  result  of  the 
introduction  of  the  integral  action  in  the  control  law. 

The  stability  requirement  is  met  whenever  the  control  parameters 

result  in  an  asymptotically  stable  system  with  adequate  relative  stability. 

To  cast  the  problem  of  the  LFC  of  a  single  area  steam  power  system  into  a 
problem  of  parameter  optimization,  the  system  dynamics  given  by  (3.1)  can 
be  written  as 


t 


f 


X  dt  +  T  AL 


(3.6) 
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Equation  (3.6)  is  the  result  of  the  direct  substitution  of  (3.2)  into  (3.1) 
Now, define  a  new  variable 


-  ci  /  V1 
to 


(3.7) 


then,  define  an  augmented  vector 


X  =  (x1,  x2,  x3,  x^) 


hence,  the  augmented  dynamic  system  can  be  written  as 


X  =  A  X  +TAL 


(3.8) 


where 


,  rT  -  (-  o,  o,  o) 


with  0^  being  3x1  row  vector  whose  elements  are  all  zeros.  The  augmented 
cost  functional  can  be  written  as 


1 - 

•tf** 

O 

-Blkl 

1 

C  3 

0 

X~QX  dt 


(3.9) 


in  which 
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qi+Ykp2 

Vi 

q2 

q3 

Vi 

q,+YK2 

Now  the  problem  is  ready  to  be  solved  by  employing  the  technique  which  has 
been  developed  in  chapter  2. 


3.2.2  The  solution  of  the  optimization  problem  of  the  LFC 

of  single  area  steam  power  system 
The  problem  of  the  LFC  of  a  SASPS  can  be  summarized  as 
follows : 


find  the  optimal  control  parameters  k_^  or 
the  cost  functional 


k  -plus-kT  which  minimize 
P  I 


J  =  \  f  XTQ(k  ,k  )Xdt 

z  '  pi 

0 

subject  to  the  dynamic  equality  constraints 


(3.10) 


X  =  AX  +  r  AL 


(3.11) 


and  to  the  inequality  constraint 


(3.12) 
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The  inequality  constraint  (3.12)  can  be  converted  into  an  equality 
constraint  [69]  by  introducing  a  new  variable 


-  x3V 


H1^P  x2  “  X3-^ 


(3.13) 


with 


x5  (tQ)=  x5 (t  )  =  0 

if  a  <  0 

if  a  >0 


in  which 


0 


H1(a)  = 


Since  the  conjugate-gradient  technique  is  employed  here  to  solve  the 
parametric  optimization  problem,  there  is  no  way  to  ensure  that  the 
constraint  will  not  be  violated  during  the  computational  procedure  of 
finding  the  optimal  parameters.  Consequently,  the  requirement  that 
*s(tf)  is  equal  to  zero  cannot  be  met.  Thus  a  penalty  is  placed  on  x^(t^)  for 
being  larger  than  zero.  Therefore  a  modified  cost  functional  can  be  written 
as  follows: 

tf 

J  =  j  I  XT  QX  dt  +  w5^(tf)  ,  (3.14) 

t0 

To  handle  a  criterion  of  the  form  (3.14),  the  system  differential 
equations  and  the  constraint  are  augmented  by  an  additional  equation  by 


introducing  a  new  variable 
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*6  =  \  XTQx  (3.15) 

with 

W  =  0 


By  defining  a  Hamiltonian 

6 

H  =  l  A,x.  (3.16) 

i=l 

and  by  applying  Pontryagin's  minimum  principle,  one  can  get  the  system 
canonic  equations  in  for  either  integral  or  proportional- 
plus-integral  control  action, 
i  -  The  optimal  integral  controller 

In  the  case  of  an  integral  controller,  the  system  dynamics 
can  be  written  as  follows 


G  ,  1  AL 

X1  M  X1  M  X3  "  M 


x 


2 


g 


(3.17) 


uo 


x6  = 


"T  I 


i=l 


2,1  n  v2 

Vi  +  2  Y(klV 


1 

2 
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with 


xi(t0)  ®  2,  ...»  6 


and  the  system  costate  equations  can  be  written  as 


qlXl  +  M  A1  +  T  X2  ”  A4 

g 


A2  " 


,  1  ,  1  ,  2  ,1 

q2X2  T  A2  ~  t  A3  “  t  j  X9 

g  t  At  t 


Tt  X3)  H1A5 


X3  " 


1  2  1  1 
qoxo  +  X  +  (—  X  -  —  x_)  H..  X^ 

Lt  J  Lt  1  Tt  J  1  ^ 


(3.18) 


*  2  1 
X.  =  -  y  x.k  +  —  k_X. 
4  4  I  T  12 

g 


with 


X.(tf)  =0,  i  =  1,  2,  3,  4. 
X5(tf)  =  2w5x5(tf) 

A6^ tf ^  =  1* 


and  the  gradient  vector  is  given  by 
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2  1 

Hk  =  ■  y  x4  kx  +  “  x4  X2  (3.19) 

1  g 

The  system  canonic  equations  have  been  programmed  for  the  computer 
with  system  parameters  given  by  [15] 

M  =  .04,  G  =  .01,  T  =  .5,  T  =  .5,  and  E  =  .03 

O  ^ 

AL=  -.005  pu 
'F  =  0.05 

and  the  weighting  coefficients  have  been  chosen  as  follows 

q.  =  1,  i  =  1,  2,  3 
Y  =  1 


The  penalty  factor,  in  the  case  that  x.(t^)  is  greater  than  zero, 
has  been  chosen  to  be 


w^  =  1000. 


Using  the  gradient  as  determined  in  (3.19)  over 


[tQ,  tf]  =  [0,  20  sec] 


conjugate  gradient-descent  has  been  accomplished  with  the  program  described 

in  chapter  2.  The  initial  guess  of  the  control  parameter  and  the 

-2 

step  size  have  been  taken  to  be  equal  to  zero  and  10  respectively. 
After  2  iterations  (13  trials) ,  the  optimal  integral  gain  parameter  has 


. 


1*1 

’ 
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been  found  to  be  as  follows 


k  =  .01703 

with  a  corresponding  cost 


J  =  0.02119 


with 


H,  I  I  =  0.16x10 


The  augmented  dynamic  system  eigenvalues  are  given  as  follows: 


-0.769,  -2.9476,  -0.2665  ±  j0.824 


Therefore,  the  optimal  integral  control  parameter  results  in  an 
asymptotically  stable  system. 


ii  -  The  optimal  proportional-plus-integral  controller 

The  problem  of  finding  the  optimal  proportional-plus-integral 

control  gain  parameters  is  considered.  Thus,  equations  number  1  and  6 
in  (3.17)  are  replaced  by 


E  1  1  .  1  , 

r  xi  -  r  x2  -  r  Vi  •  r  V* 

g  g  g  g 


and 
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=  1  I 

i=l 


2,1  „ 

qi  i  +  2  Y(klx4 


+  k  x.) 
P  1 


to  obtain  the  system  dynamic  equations  in  the  case  of  3 
proportional-plus-integral  control  action.  Equations  1  and  A  in 
(3.18)  are  also  replaced  by 


xi  -  - 


qixi 


-  Yk  x  -  yk  k  x. 
pi  p  I  A 


A1  + 


E_ 

T 


g 


X2  ' 


and 


X4  =  -  YkIx 4  -  Ykpklxl  +  T  kix2 

g 

to  get  the  system  costate  equations. 

In  this  case  the  gradient  vector  components  are  given  by 


Hk  =  -  Y(kpXl  +  kIxlx4)  +  T  x1A2 
P  g 


and 


-  y(klX2  +  k  xxx4)  +  i-  x4A2 

g 


Using  the  gradient  as  determined  in  (3.19)  over  the  same  time  interval 
with  the  same  system's  parameter  and  load  disturbance,  the  conjugate-gradient 
descent  has  been  accomplished  with  the  algorithm  described  in  chapter  2.  The 
initial  guess  of  the  control  parameters  k  ,  k^  and  the  step  size  have  been 

_3 

taken  simply  to  be  zero,  zero,  and  10  .  After  8  iterations  (23  trials), 


•  ; 
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the  optimal  control  parameters  are  as  follows 

k  =  0.0859 
P 

k  =  0.02177 

and  the  corresponding  cost 

J  =  0.696xl0-2 


with 


H.  =  0.2926x10 

ii  k  i  i 

The  eigenvalues  of  the  augmented  system  are  given  as  follows 


-3.769,  -0.1484  ±  jl.7641,  -0.18368 


Therefore,  the  optimal  proportional-plus-integral  controller  results 
in  an  asymptot ically  stable  system.  However,  in  designing  a  control  system, 
it  is  not  only  required  that  the  system  be  stable  but  also  it  is  necessary 
that  the  system  has  adequate  relative  stability.  Here,  the  optimal 
system  relative  stability  is  unsatisfactory  since  it  has  a  pair  of  dominant 
complex  -conjugate  closed  loop  poles  near  the  jw  axis.  However,  the 
problem  of  designing  an  optimal  output  feedback  controller  with  satisfactory 
relative  stability  can  be  achieved  by  modifying  the  system  structure.  This 
can  be  realized  by  addition  of  a  zero  to  the  open-loop  transfer  function 


' 
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which  has  the  effect  of  pulling  the  root  locus  to  the  left.  Physically, 
the  addition  of  a  zero  in  the  feedforward  transfer  function  means  the 
addition  of  derivative  control  to  the  system.  Therefore,  the  case  of 
studying  the  addition  of  derivative  to  the  proportional-plus-integral 
controller  is  considered  next. 


“The  optimal  proport ional— plus-integral— plus— derivat ive  controller 
The  control  law  in  this  case  can  be  written  as  follows 

tf 

u  =  -  Vi '  ki  1  xidt  -  V-i 

to 


-I  V 


M  kdX3 


kIX4  +  j?  kd  &L 


(3.20) 


By  applying  the  procedure  adopted  in  section  (i)  and  (ii),  the  system 
dynamics  equations  are  the  same  as  in  eqn.  (3.17)  except  x?  which  is  replaced 
by 


G_ 

MT 


g 


V  xi 


MT 


g 


X3  -  T~  klX4 


+ 


MT 


ALk 


g 


Also  x^  is  replaced  by 


=  l  q  .x2  +  [(k  -  ^  k,)x  +  ^  -  k  x  +  ^  k  AL] 

6  ^  ^ii  1  p  Md  l  M  d  3  14  Md 


The  system  costate  equations  are  also  the  same  as  in  eqn.  (3.18)  except 
is  replaced  by 


' 
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II 

tH 

•  r< 

*  qixrYt(kP  - 1  kd)xi  +  5  Vs  +  kix4  -  5  VL1(kP  -  jj  kd) 

+  m  Ai  +  (f  +  r  kP  -  sr  kd>  ^  \ 

Also  A^  is  replaced  by 


q3X3  “  M  kpt(kp  “  M  kd)xl  +  M  kdX3  +  kIX4  "  M  kdAL-* 

M  X1  +  MT  kdX2  +  7~  X3  +  T~  (T~  X2  "  f  x3)HlX5 
g  t  t  t  At 

Moreover  A^  is  replaced  by 


v  • 

i'¬ 

ll 

-  r  k  [(k  -  -jj  k  )x  +  ^  k  x  +  k  x.  -  k  AL] 

I  p  Mdl  MdJ  14  M  d 

-  k  A 

T  I  2 

g 

The  gradient  vector  components  are  given  by 


P 

-  rxlI(kp  -  1  kd)xl  ±  5  kdx3  +  kIx4  -  i  kdiLl 

+  FX1X2 

g 

H.  = 
kt 

X 

-  rx. [(k  -  ^  k  )x  +  ^  k  x  +  k  x,  •-  ^  k  AL]  +  x  A. 

4  p  Mdl  M  d  3  14  Md  T  42 

g 
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H 


k 


d 


1  AL 

M  X3  ‘  M 


(3.21) 


Using  the  gradient  given  by  equations  (3.21)  over  the  time  interval 
U0,tf]  =  [0,20sec] 

the  conjugate  gradient  is  found  with  the  technique  described 
in  chapter  2.  The  initial  guess  of  the  control  parameter  vector  and  the 
step  size  have  been  simply  chosen  equal  to  0.107,  0.077,  0.08276  and  10“ 
respectively.  After  30  iterations  (62  trials),  the  optimal  control 
parameters  have  been  found  as  follows 

k  =  0.111 
P 

k  =  0.0790 

k,  =  0.0863 
d 

and  the  optimum  cost  is 

J  =  0.18355xl0"2 


with 


H. 


K 


0.0102 


i  'i 
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The  system  s  eigenvalues  corresponding  to  the  optimal  control  vector 
are  given  as  follows 

-0.724  ±  j0.6226,  and  -1.4  ±  J2.588 

Therefore,  by  employing  a  proport ional-plus-integral-plus— derivative 
controller,  a  system  with  adequate  relative  stability  as  well  as  less 
state  deviations  has  been  established. 

Quite  often  in  designing  a  controller  for  a  higher  order  system  the 
control  parameters  are  adjusted  so  that  there  will  exist  a  pair  of 
dominant  complex— conjugate  closed-loop  poles.  The  presence  of  such  poles  in 
a  stable  system  reduces  the  effect  of  such  nonlinearities  as  dead  zone 
and  backlash.  Therefore,  a  proportional-plus-integral-plus- 
derivative  controller  will  be  designed  next  so  that  the  system  will  have  a 
preassigned  dominant  pair  of  complex -conjugate  poles  and  the  rest  of  the 
poles  will  be  preassigned  negative  real. 

iv  -  Optimal  proportional-plus-integral-plus-derivative  controller  by 

pole-assignment 

The  characteristic  equation  of  a  single  steam  area  power  system  can 
be  written  as  follows 

5  4  . 

I  a  s  _1  =  0  '  (3.22) 

i=l  1 


where 
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a3  MT  +  MT  +  T  T  +  MT  T  kd 
8  t  t  g  t  g 


=  G  +  E  1 _ 

a4  MT  T  MT  T  p 
g  t  g  t  F 


a5  MT  T  kl 
g  t 


(3.23) 


B>  substituting  the  system  parameters  in  (3.23),  one  can  write  (3.24) 
as  follows 

s  4  +  4.25s3  +  (5  +  100kd)s2  +  (4+100kp)  s  +  lOOkj  =  0  (3.24) 

The  problem  posed  is  to  find  the  optimal  control  parameters  k  ,  k 

P  I 

an^  kd  reclu^re^  to  position  the  four  poles  of  the  plant  plus  the  integrator 
at 


“  *5  ±  jO ,  -  1.25,  and  -  2  f  <j>  >  o 

By  formulating  the  desired  characteristic  equation  and  comparing  it 
to  the  closed-loop  one,  one  can  get 

k  =  -  0.039375  +  3.25kJ 
P  d 

(3.25) 


k_  =  -  0.01875  +  2.5k, 
I  a 
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To  find  the  optimal  pole  assignment  control  parameters,  one  can 
adopt  the  same  procedure  as  in  section  (i) .  This  results  in  a  system 

canonic  equations  similar  to  the  one  given  in  section  (iii)  except  that 
is  replaced  by 


*2  “  -  [f"  +  Y"(-  0-039375  +  3.25kd)  -  *L-  kdl* 

g  g  g 

-  T"  x2  -  ST  x3  -  T  (-  °'01875  +  2-5V  x4  +  i 

g  g  g 


d'  "4  '  MT  ALkd 
g 


Also  A^  is  replaced  by 


X1  =  -  -  Y[(-  0.039375  +  3.25kd  -  |  kfl)  ^  ±  k^ 

+  (-  0.01875  +  2.5  kd)x4  ^  kd  AL][3.25  ^)kd  -  0.039375] 


+  f  *i  +  if-  +  T~  (~  °*039375  +  3.25kd) 
g  g 


MT  kd  ^  X2  " 


^  is  replaced  by 


*3  =  -  93x3  -  5  0.039375  +  3.25kd  -  £  V  ^  +  £  k  x3 


+  (-  0.01875  +  2.5k,)x,  -  ^  k0  AL] [-  0.039375  +  3.25k,] 

d  4  M  d  d 


1,  J  ,  ,  1  ,  ,  2  A 

"  M  A1  MT  kd  A2  +  T  A3  +  T  (T  X2 

g  t  t  t 


Tt  X3)  H1A5 


Finally  A^  is  replaced  by 


. 

'■ 
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-  -  V4  -Tic-  0.039375  +  3.25^  -  £  k,)  x,  4  ±  krf  x3 
+  (-  0.01875  +  2.5kd)x4  -  i  kjAL] [-  0.01875  +  2.5kj] 
-  Y~  (-  0.01875  +  2.5kd)X2 


The  gradient  with  respect  to  kd  is  given  as  follows 


H. 


-  Y  [(-0.039375  +  3.25k,  -  £  k(J)  x,  +  ±  kd  x3 


+  (-°-°1873  +  2'5V  x4  -  5  VL1  [<3.25  -  §)  x,  +  ±  x3 


+  2.5  x  -  ~  AL]  +  [-^1  x  +  —  +  —  X  -  — 1A 
4  M  J  1  T  A1  MT  T  x/.  j  A 


T  1  MT  T  4  MT  J  2 
8  g  g  g 


(3.26) 


By  choosing  the  initial  guess  of  the  control  parameter  k,  and  the  initial 

d 

-2  -5 

step  size  equal  to  1x10  “  and  10  respectively,  the  optimal  control 
parameter  is  found  to  be  as  follows 

kd  =  0.096 

and  the  optimal  cost 

J  =  0.147  x  10~2 
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and  the  gradient 

|  |H  ||  =  0.67  x  1C)"2 
kd 

In  conclusion,  the  problem  of  designing  an  optimal  or  subopt imal 

controller  for  the  problem  of  LFC  of  single  area  steam  power  system  has 

been  considered.  In  the  case  of  employing  integral  and/or  proportional 

control  action,  the  optimal  control  parameters  k  and/or  k  have  been 

1  P 

found  to  be  independent  of  the  load  changing.  In  the  case  of  adopting 
a  PID  control  strategy,  the  optimal  control  parameters  have  been  found 
to  be  strongly  dependent  on  the  load  changing.  The  system  response  for 
different  controllers  in  the  case  of  disturbing  the  area  by  -  0.005  p.u., 
is  shown  in  fig.  [2].  Fig.  [3]  shows  the  system  response  for  the  optimal 
PID  controller  and  the  suboptimal  one.  They  emphasize  the  figure  of  merit 
of  adopting  modern  control  theory;  rather  than  the  classical  one,  on 
designing  system  controllers. 

In  spite  of  the  fact  that  the  optimal  or  the  suboptimal  PID  controllers 
give  better  performance  than  the  optimal  integral 

and/or  proportional  controllers,  it  is  recommended  that  the  optimal  PI  controller 
be  used  for  the  following  reasons: 

-  The  optimal  PI  control  parameters  are  independent  of  the 
load  changing 

-  It  is  simpler  to  design  the  hardware  of  the  PI  controller 
rather  than  that  of  the  PID  one 

-  The  PID  controller  might  produce  an  undesirable  system  performance 
in  the  case  that  the  area  control  signal  experiences  a  noise 


. 

' 
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component  in  its  measurement. 

3.3  Load  frequency  control  of  single  area  hydro  power  system 
3.3.1  The  problem  formulation 

A  typical  linearized  model  of  a  hydro  area  power  system;  (SAHPS)  is 
shown  in  fig.  (4)  [15] .  The  system  dynamics  can  be  put  in  the  state 
space  form  as  follows: 

X1  =  W  +  Blu  +  riAl  (3.27) 


where 


(Aw,  Ax  ,  A  ,  Ap  ) 
g  g  g 


in  which 


A 

g 


the  deviation  in  the  gate  position  in  pu  power 


and 


C*  =  (1,  0,  0  ,  0)  ,  B*  =  (1.  o,  0,  0) 

rT  =  ^  o,  o) 
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Fig.  (2-a)  Frequency  deviation-time 


Fig.  ( 2— b )  Output  power  deviation-time 
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o 


Fig.  (3-a)  Frequency  deviation-time 
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Fig.  (3— b )  Output  power  deviation-time 


Comparison  between  the  optimal  and  the 
suboptimal  P-I-D  controller  on  the  dynamic 
response  of  single  area  power  system 
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The  control  law  is  specified  to  be  of  integral  or  proportional— 
plus-integral  action.  Therefore,  the  control  law  can  be  written 
in  general,  in  terms  of  the  system  states  as  follows: 

tf 

u  =  "  kp  X1  -  kI  I  ^dt  (3.28) 

t0 

To  cast  the  problem  in  the  normal  state  space  form,  a  new  variable 
is  defined  as  follows: 


then,  define  an  augmented  vector  X  by 


(3.29) 


hence,  the  augmented  dynamic  system  can  be  written  as  follows: 


X  =  AX  +  T  AL 


(3.30) 


1 
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where 


and 


rT  =  (-  i  ,  o,  o,  o,  o) 


The  problem  posed  is  to  find  the  optimal  values  of  k  and  (k  k  ) 

I  pi 


which  minimize  the  cost  functional 


J  =  j  /  (XTQX  +  ru2)dt 


(3.31) 


such  that  the  dynamic  system  given  by  (3.30)  is  asymptotically  stable 


3.3.2  The  optimal 
To  handle  a  criterion 
equations  are  augmented  by 
defined  by 


solution 

of  the  form  (3.31),  the  system  differential 
an  additional  equation.  This  equation  is 


■ 
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(3.32) 


with 


W  =  0 


By  defining  a  Hamiltonian 


6 


(3.33) 


then,  by  applying  the  Pontryagin  minimum  principle,  one  can  get  the 
system  canonic  equations.  By  solving  these  equations  as  described  in 
chapter  2,  one  can  get  the  optimal  control  parameters. 

i  -  The  optimal  integral  controller 

In  the  case  of  employing  an  integral  controller,  the  control  law 
which  represents  the  speed  changer  commands  can  be  written  as  follows: 

u  =  -  \  I  x. jdt  (3.34) 

*0 

Therefore,  following  a  step  load  change,  a  new  steady  state  is  attained 
after  the  speed  changer  command;  u,  has  reached  constant  value.  But  this 
evidently  requires  that  the  integrand  in  (3.34)  be  zero,  i.e., 


(3.35) 


which  is  the  basic  proposition  of  the  LFC  in  single  area  systems.  The 
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application  of  the  Pontryagin  minimum  principle  results  in  the  following 
system  canonic  equations 


•  _  G  1  AL 

X1  ~  "  M  X1  M  X4  ~  M~ 


.  _  E  1  1  , 

X2  T  X1  “  t  X2  ~  T  ^1  X5 

8  g  g 


.  _  1  1 

X3  Tt  X2  ~  T  X3 


x4  T  X2  +  (D  +  T  }  X3  “  D  x4 


(3.36) 


xc  =  x 


1  r  2  1 

~2  \  qiXi  +  2  r(kIx5) 
1=1 


with 


Xi(t0)  =  °  i  =  1,  2, 


6 


and 


X 

1 


Vi  + 1  h 


X2  "  - 


q2x2 


g 


A,  -  — 

2  Tt 


\ 

3  xt  4 


q3X3 


57 


s  h  + 1  \ 


.2 

r  kIx5 


X 


6 


0. 


(3.37) 


with 


A1(tf)  =0  i  =  1,  2,  5 

VV  ■ 1 


with  a  gradient 


2. 

rx5kI 


(3.38) 


The  system  canonic  equations  have  been  programmed  to  be  solved  with 
system  parameters  given  by  [15] 


M  =  0.03,  G  =  0.003,  T  =  1,2,  T 

g  t 


0.5,  D  =  0.5,  and  E  =  0.013 


and  a  step  load  disturbance  given  as  follows 


AL  =  4-  0.005pu 


and  the  weighting  coefficients  have  been  chosen  as  follows 
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<3^  1,2,  . .  . ,  4 

r  =  1 

Using  the  gradient  as  in  (3.38)  over 
[tQ,tf]  =  [0,20sec] 

with  an  initial  guess  for  the  control  parameter  and  the  step  size  equal 

-4 

to  zero  and  10  respectively,  the  conjugate-gradient-descent  has  been 
accomplished  in  2  iterations  (9  trials).  The  optimal  integral  gain 
parameter  has  been  found  to  be  as  follows 

k  =  0.00355 

and  the  corresponding  cost  is 

J  =  0.19211 

with  a  gradient  norm 

I |H,  ||=  0.242  x  10~3 

V 

The  augmented  dynamic  system  eigenvalues  are  given  as  follows 

-  0.1586  ±  j0.494S  and  -  3.253  ±  j0.1444 
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Therefore,  the  optimal  integral  controller  results  in  an  asymptotically 
stable  system. 


ii  ~  The  optimal  proport ional— plus-integral— controller 

In  the  case  of  adopting  a  proportional-plus-integral  control  action, 
the  control  law  can  be  written  as  follows 

u  =  "  kDxi  -  kT  /  xidt  (3.39) 

F  t 

0 

Following  a  step  load  change,  at  steady  state,  the  speed  changer  command 
is  constant.  This  requires  that  the  integrand,  the  frequency  deviation, 
in  (3.39)  be  zero. 

The  augmented  system  dynamic  equations  is  the  same  as  in  section 
(i);  eqn .  (3.36),  except  that  x^  is  replaced  by 


k  ) 
P 


1  1  i 

—  x2  -  — 

8  g 


and 


is  replaced  by 


=  t  I 


i=l 


qixi  +  \  r(kpxl  +  kIx5)2 


with  the  same  initial  conditions. 


Now,  by  defining  a  Hamiltonian  as  follows 


H-  I 


i=l 


w 


. 
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then,  by  applying  the  Pontryagin  minimum  principle,  one  can  get  the 
costate  equations  as  in  section  (i);  eqn.  (3.37),  except  that  X^  is 
replaced  by 


xi "  - 


qlXl  '  rkp< kpxl  +  kIX5)  +  M  X1  + 


+  T  kp^  X2 
g  g 


and  the  gradient  vector  components  are  as  follows 


Hk  =  -  rxl(kpxl  +  W  +  T  X1  x2 
p  g 


and 


H.  =  -  rxc(k  x-  +  k  x,.)  +  —  x  X 
k.  5  p  1  15  T52 

1  g 


(3.40) 


Using  the  gradient  given  by  equations. (3 .40)  over  the  time  interval 


[tQ,  t  ]  =  [0,  20sec ] 

The  control  parameters  were  found  in  9  iterations  (30  trials) .  The 
initial  guess  of  the  control  parameter  vector  and  the  step  size  were  taken 
equal  to  zero  and  .  1  x  10  respectively.  The  optimal  control  parameters 
are  given  as  follows 


k  =  0.0124 
P 

k  =  0.00347 

and  the  corresponding  cost  is  given  as  follows 
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J  =  0.1355 

with  a  gradient  norm 

| |H  | |=0.2976  x  10"A 

The  augmented  dynamic  system  eigenvalues  are  given  as  follows 

-  0.0853  ±  j 0 .714  and  -  3.405  ±  j0.843 

Therefore,  the  optimal  proport ional— plus— integral  controller 
results  in  an  asymptotically  stable  system.  The  optimal  proportional- 
plus-integral  controller  (0P1C)  is  better  than  the  optimal  integral 
controller  (OIC)  in  the  sense  that  the  system  overall  cost  is  smaller.  On 
the  other  hand,  the  system  relative  stability,  in  the  case  of  employing  the 
OIC,  is  better  than  the  system  relative  stability  in  the  case  of 
employing  the  OPIC.  Both  the  01  and  OPI  control  parameters  have  been 
found  to  be  independent  of  the  load  changing.  Fig.  (5)  shows  the  system 
response  in  the  case  of  a  -0,005  p.u.  load  change. 


o 


Fig.  (5-a)  Frequency  deviation-time 


0/P  POWER  DEV. 

0*0076  -0.0C6R  -0.00S8  -0.0016  0.0006 
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Fig 


(5-b)  Output  power  deviation-time 
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I  controller 


P-I  controller 


Fig.  (5-c)  Area  control-time 


Dynamic  response  of  a  single  hydro  area  power  system 


Chapter  4 

Load  Frequency  Control  of  Multi  Area 
Interconnected  Power  System 


4.1  Introduction 

A  multi  area  interconnected  power  system  is  one  in  which  each  area 
of  the  system  is  expected  to  adjust  its  own  generation  to  absorb  its 
own  load  changes.  There  is  no  attempt  to  have  a  complete  coordination 
between  all  the  areas  of  the  interconnected  system.  Each  area  deals 
primarily  with  its  adjacent  neighbors.  Tie-line  power  flow  between 
adjacent  areas  are  scheduled  and  maintained. 

The  load  frequency  control  of  an  area  of  an  interconnected  power  system 
relies  on  an  operating  schedule  which  is  often  prepared  24  hours  ahead 
of  time.  This  schedule  indicates  the  expected  demand  profile  of  the  area 
as  well  as  the  area  commitment  to  its  adjacent  areas.  The  area  commitment 
is  the  tie-line  power  interchange  which  should  be  maintained  at  a  certain 
point  of  time.  These  values  are  fed  to  the  area  controller  as  governor 
set  points. 

In  the  normal  mode  of  operation,  the  area  prevailing  frequency  and 
tie-line  power  are  compared  to  the  area  set  point,  scheduled  frequency  and 
tie-line  power,  every  2-5  seconds.  An  area  control  error  (ACE)  is  computed 
according  to  the  following  equation 

ACE  =  AP  +  BAf 


where 

AP  =  area  net  interchange  tie-line  power  deviation 
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Af  =  frequency  deviation 
3  =  area  bias-setting 


The  ACE  is  converted  to  a  train  of  pulses  of  raise/lower  commands 
which  drives  the  speed  changer  of  the  governing  system  in  accordance  with 
the  sign  of  the  ACE.  A  negative  ACE  causes  generation  to  increase  and 
vice  versa.  The  speed  governor  shaft  position  will  change  in  steps.  At 
the  end  of  the  control  period  the  speed  governing  characteristic  will  be 
shifted  to  a  new  position,  a  new  set  point,  which  results  in  matching  the 
load  demand  with  the  load  consumption. 

Therefore,  the  governor  set  points  are  continuously  altered  to  keep 
area  frequency  and  tie-line  power  at  their  scheduled  values.  The  transition 
from  one  set  point  to  another  represents  a  step  change  in  the  load  demand 
which  itself  introduces  a  frequency  transient  and  tie-line  power  disturbances 
before  settling  down  to  their  scheduled  values.  Therefore,  it  is  necessary 
to  design  the  area’s  controllers  with  control  gain  parameters  such  that  the 
area’s  transients  are  minimized  as  well  as  the  area’s  steady  state 
specifications  being  met. 

The  problem  of  designing  area  frequency  regulators  using  conventional 
and  modern  control  theory  has  been  the  subject  of  numerous  studies  since 
the  earliest  days  of  generation  [5-30].  A  significant  group  of  results 
covering  various  aspects  of  this  problem  are  available  in  the  literature, 
(some  of  the  relevant  references  are  given  at  the  end  of  this  dissertation) . 
The  conventional  LFC  approach  often  employs  the  tie-line  bias  concept 
in  designing  the  area's  controller  which  has  proportional  and/or  integral 
action.  This  type  of  control  is  extensively  used  in  practice,  inspite  of 
all  techniques  which  have  been  developed,  recently,  employing  modern 
control  theory  [5].  The  reason  for  that  is  the  fact  that  most  of  these 
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recent  techniques  have  developed  linear  feedback  controls  which  are 
functions  of  all  the  system  states  as  well  as  system  disturbances  [6],  [7]. 
Therefore,  it  was  necessary  to  design  an  observer  to  realize  these  kinds 
of  controls  [8],  Once  an  observer  is  introduced  into  the  system, the  cost 
will  be  increased  and  the  control  is  no  longer  optimal  [9] .  Another 
important  reason  is  that  a  control  which  depends  upon  all  the  system  states 
needs  that  some  of  these  states  to  be  telemetered  since  the  areas  of  inter¬ 
connected  power  systems  are  usually  spread  over  large  geographical 
territories.  This  is  why  in  practice,  control  engineers  prefer  to  use 
the  conventional  control  to  the  advanced  one  inspite  of  the  fact  that  the 
latter  improves  system  transient  performance. 

Calovic  [10] ,  [11]  was  the  first  to  narrow  the  gap  between  the 
conventional  and  the  advanced  LFC  by  adopting  Porter’s  idea  [12]  of 
introducing  the  integral  of  the  area  output  states  as  a  part  of  the  control 
design.  However,  the  proportional  part  of  the  control  was  a  function  of 
all  system  states  and  an  observer  was  needed  to  realize  such  a  control. 
Here,  modern  control  theory  is  implemented  to  find  the  optimal  control 
parameters,  of  the  integral  and  the  proportional  parts  of  the  control  law 
in  a  systematic  way  in  order  to  meet  the  systems  transient  and  steady 
specifications,  as  well  as  to  conform  to  what  control  engineers  practice  in 
the  real  world  of  load-generation  control. 

4.2  Decentralized  optimum  load  frequency  control  of  interconnected 

power  systems 

The  case  considered  in  this  dissertation  is  that  of  two  mixed  area 
interconnected  power  systems  (TMAIPS)  which  may  be  regarded  as  a 
representation  of  a  particular  area  which  undergoes  a  disturbances,  connected, 
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via  a  tie-line,  to  the  rest  of  the  system.  The  disturbed  area  is  assumed 
to  be  subjected  to  a  step  load  change  AL.  The  justification  of  assuming 
that  interconnected  power  systems  are  subjected  to  deterministic  load 
disturbances,  inspite  of  the  fact  that  they  are  normally  subjected  to 
random  disturbances,  is  that  modern  LFC  systems  are  designed  with  filters 
[13]  in  order  to  remove  the  purely  random  portion  of  the  regulating 
responsibility  which  generation  can  not  follow,  leaving  the  deterministic 
components  for  the  system  units  to  follow. 

Once  an  area  has  experienced  a  load  change;  e.g.  load  increase,  a  chain 
of  events  ensues  as  follows: 

( i)  The  disturbed  area  speed  begins  to  decrease  as  the 
increased  load  demand  is  supplied  from  the  stored  kinetic 
energy 

(ii)  The  tie-line  phase  angle  increases  and  unscheduled  tie-line 
power  flows  to  the  disturbed  area 

(iii) The  remote  area  speed  begins  to  decrease,  as  the  disturbed 
area  unscheduled  tie-line  power  represents  a  remote  area  load 
increase 

(iv)  The  governors  on  both  areas  sense  the  change  in  the 
speed  and  each  area  responds  in  proportion  to  its  natural 
governing  characteristic 

(v)  The  supplementary  regulators  come  into  action  either 
after  the  governors  have  ceased  operating  or  during 
their  operation. 

The  approach  which  will  be  adopted  in  designing  the  supplementary 
regulators  depends  on  whether  they  come  into  action  before  or  after  the 
governors  have  stopped  operating.  Load  frequency  control  operation  in 
which  the  governors  cease  operating  before  the  action  of  the  supplementary 
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regulators  will  be  referred  to  as  non-interactive  governor-supplementary 
regulator  load  frequency  control .  The  case  in  which  the  supplementary 
regulators  come  into  action  during  the  operation  of  the  governors  will 
referred  to  be  as  interactive  governor-supplementary  regulator  load 
frequency  control. 

Mathematical  Model 

A  linearized  model  of  two  mixed  area  interconnected  power  system 
(TMAIPS)  may  be  written  as  [8]: 


*1  1 

X  =  A^X  +  B1U  +  riAL 
Y  =  C^X1 


(4.1) 


in  which 


i  1  i  '  !  Iii 

-  TJ-'  0  0  I  0  j  0  10'  0*0 

”l I  I  I  I  I  I 


where 


the  frequency  deviation  of  area  1 
the  frequency  deviation  of  area  2 
the  power  angle  deviation 

the  deviation  in  the  governor  position  of  area  1 

the  deviation  in  the  governor  position  of  area  2 

the  deviation  in  mechanical  output  of  area  1 

the  deviation  in  the  gate  position  of  area  2 


x2  =  the  deviation  in  mechanical  output  of  area  2 


4 
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^•2*1  The  optimal  non-interactive  governor-supplementary  regulator 
In  this  mode  of  LFC  operation  the  interconnected  system  settles 
down  at  a  new  frequency  common  to  the  areas  of  the  interconnected  system  and 
a  new  tie-line  flow  before  the  supplementary  regulator  comes  into  action 
the  case  of  a  load  increase;  e.g.  in  area  1,  AL,  the  system  common  frequency 
is  given  [14] 


f  = 


fS  + 


-AL 


2t:[(G1+E1)  +  (G2+E2)] 


and  the  tie-line  power 


-AL(  G  +E  ) 

p  =  ps  +  - - - _ - 

t12  t12  (G1+E1)+(G2+E2) 


(4.2a) 


(4.2b) 


in  which 


f  =  the  scheduled  frequency 
pS 

t^2  =  the  scheduled  tie-line  power 


It  is  desired  to  design  an  optimal  supplementary  regulator  with 
proportional-plus-integral  control  of  the  area  control  error  such  that 
system  frequency  and  tie-line  power  are  maintained  at  their  scheduled  values. 
In  terms  of  the  system’s  state  variables  the  control  law  can  be  written  as: 


u(t)  =  -  K  ace(t)  -  K  /  act(t)dt 
”  1  t 

0 


(4.3) 


in  which 


*. 
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ace(t)  =  D  x^t) 


where 


S  =  dia§  (kpl’V 

Kj  =  diag  (ku,kI2) 


0 

1 

0 

0 

0 

0 

0 

0 

B2 

-1 

0 

0 

0 

0 

0 

The  control  parameters  and  can  be  expressed  mathematically  as: 


k  .  =  0  i  *  1,2 

pi 

(4.4) 

l  =  0  i  -  1,2 


The  problem  posed  is  to  find  the  optimal  control  parameters  k 
and  k  ,  i  =  1,2,  which  minimize  the  system’s  transient  and  the  control 
action  such  that  the  system’s  steady  state,  dynamic  limit  and  area’s  de¬ 
centralization  are  met. 

Minimization  of  the  system’s  transient  and  the  control  effort  can  be 
accomplished  by  minimizing  the  cost  functional 


1 

2 


I  f  (X] 


Qlx 


+ 


UTRU)dt 


(4.5) 
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with  respect  to  the  control  parameters  and  subject  to  the  dynamic 
constraints  given  by: 


(4.6) 


in  which  ,  B^,  T  are  defined  in  (4.1). 

Here  and  R  are  diagonal  positive  definite  matrices  of  dimension 
8  and  2  respectively. 

The  dynamic  requirement  is  met  when  the  steam  area  generation  limit 
is  kept  within  \p  p.u.  mw/sec .  the  generation  rate  limit  is  not  fixed  and 
varies  from  unit  to  unit  (60  ip  =  2-5) .  It  can  be  expressed  mathematically 
in  terms  of  the  system’s  state  variables  as: 


(4.7) 


in  which 


1 


1 


0  0) 


0 


F1  =  (  0  0  0 


T 


T 


tl 


tl 


Ideally,  to  ensure  the  nonintervention  principle  of  the  inter¬ 
connected  load  frequency  control  operation,  that  is,  power  generation  must 
be  altered  only  in  the  disturbed  area;  area  frequency  bias  setting  equal 
to  the  area  natural  governing  characteristic  must  be  employed  in  the  area 
control  error  equation.  When  this  condition  is  fulfilled,  the  control 


error,  in  the  case  of  employing  non-interactive  governor-supplementary 
regulator,  of  the  remote  area  will  equal  zero  and  its  supplementary 
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regulator  remains  inoperative. 

To  find  the  optimal  control  parameters  and  K^,  the  system’s  dynamics 
given  by  equation  (4.6)  can  be  written  as: 


X1  =  ^X1  -  B^DX1  -  BpKpD  /  x1dt  +  I^AL  (4.8) 

to 

Equation  (4.8)  is  the  result  of  the  direct  substitution  of  equation  (4.3) 
into  equation  (4.8)  .  Now  define  a  new  variable  z  by 

tf  , 

z  -  D  /  X  dt  (4.9) 

t0 

then,  define  an  augmented  vector  X  such  that 


x .  = 

l 


1 
x . 

l 


j  =  1,2 


(A. 10) 


hence,  the  augmented  dynamic  system  and  the  cost  functional  can  be  written 
as 


X  =  AX  +  F  AL 


(A. 11) 


J 


(A. 12) 


where 


A  = 


•» 
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with  _0  being  the  2x2  matrix  whose  elements  are  all  zeros, 


T  T 

r  =  (r*,o,o) 


and 


Q  = 


- L _ _ 

K^RKpD  I  Q2+KpRKI 


in  which  is  a  diagonal  positive  definite  matrix  with  dimension  2. 

The  inequality  constraints  given  by  equation  (4.7)  can  be  written  as 


FX  <  \p 


(4.13) 


in  which 

FT=  (F*  0  0) 

The  optimal  solution 

The  inequality  constraints  given  by  equation  (4.13)  can  be  converted 
into  equality  constraints,  as  in  Chapter  3,  by  defining  a  new  variable 
x^  such  that 

0  FX  _<  0 

x  =  (FX) ^  H(FX),  H(FX)  =  {  (4.14) 

11  1  FX  >  0 


To  handle  a  criterion  of  the  form  (4.12),  the  system  differential 
equations  and  the  constraint  are  augmented  by  an  additional  equation  as 
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follows 


1  T 

x12  =  J  X  QX 


(A. 15) 


Define  the  Hamiltonian 


12 


H-  I 


i=l 


X  .x . 

l  l 


(A. 16) 


then,  by  applying  Pontryagin's  minimum  principle,  one  can  get  the  syst 
canonic  equations.  In  the  case  of  employing  proportional-plus-integral 
control  action,  the  system  dynamics  can  be  written  as  follows 


em 


X1  = 


fl  1_  1  AL1 

^  X1  M  X3  M  X6  "  M 


.  _  ^2  1  1  _  AL2 

X2  M2  X2  M2  X3  M2  X8  M2 


X3  T12  X1  T12  X2 


x4  = 


E1  1  3lkPl  kPl  kIl 

— -  X.  -  — —  x.  -  — — -  X.  -  — -  X-  — —  x_ 

T  ,  1  T  ,  A  T  ,  1  T  ,  3  T  ,  9 


gl  gl 


gl 


gl  gl 


x5  = 


2  1  B2kp2  ,  kP2  kI2 

x„  -  — —  X _  -  — -  x„  4-  — —  X^  -  — x 


T  '2  T  0  5  T  0 

g2  g2  g2 


2  T  0  3  T  9  10 

g2  g2 


1  1 

X ,  =  — —  x.  -  — —  X, 
6  Ttl  4  Ttl  6 


X-,  = 


7  T 


t2 


X5  T 


t2 


X_  =  — 


X.  + 


Tt2  "5  (Tt2+D2)  7 


x_,  - 
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Xg  =  B^x^^  +  x^ 


x10  62X2  X3 


X11  (T  XA  '  T  X6'W  H1(T  X4  '  T  . 

tl  tl  tl  tl 


x12  =  2  Vi  +  I  rl[kPl(BlXl  +  X3>  +  kl1  9 

1=1 


xj 


+  2  r2[kP2(62X2  "  X3}  +  kl2X10^ 


(4.16) 


with 


xl(t0)  =  x2<t0)  = 


AL 


(e1+g1)+(e2+g9) 


x3(t0) 


al(e2+g2) 


(e1+g1)+(e2+g2) 


W  =  X6(t0)  =  -  ElXl(tO} 


x5(t0)  x7^tQ')  x8^t0^>  ~E2X2^t0^> 


xi(t0)  =0  i  =  9,10, . . ,12 


and  the  system  costate  equation  can  be  written  as 


1  =  -qixrrl6lkpl(3lkPlXl  +  kPlX3  +  kUX9) 


G  El"*"^lkpl 

+  S7h-T12A3  +  (-T— >  X9 

1  g1 


78 


*2  '  *  q2x2  -  r262kP2(62kP2X2  "  kP2x3  +  kI2x10) 


G2  E2  +  32kP2 


M2  2  12  3 


g2 


X3  =  -  q3X3  -  rlkPl(BlkPlxl  +  kplx3  +  knx9) 


+  r2kp2<'S2kp2x2  kP2x3  +  kI2x10)  +  M.  X1 


kpi  k. 


1  "PI  "P2 

mT  A2  +  F7  X4  "  77  X5  "  x9  +  A10 
2  g1  g2 


^4  ~  ~  ^4^  +  x —  ^4  “  t —  ^6  ~  T —  —  X4  “  x —  X6~^  Hl^ll 


gl  tl  tl  tl 


tl 


;  _  ,  1  1  %  ,  2 
a,  =  -  qcx  +  — —  A  -  — —  A_  +  — —  A0 
5  5  5  T  _  5  T  „  7  T  „  8 


g2  t2 


t2 


x6  =  -  Ve  -  t  xi  +  7  x6  +  T7  (7T  x4  -  T~  V*}  H1X11 

1  tl  tl  tl  tl 


•  1  ?  2 

A  =  -  q_x^  +  ^ —  A  -  7 —  +  ^-)  A0 

7  7  7  Tt2  7  Tt2  D2  8 


X8 


q8x8  ‘  M2  X2  +  D2  X8 


x9  =  -  q9x9  -  r1kn(61kplx1  +  kplx3  +  knx9)  +  ^  X4 

g1 


X10  “  q10X10  r2kl2(62kp2X2  kp2X3  +  kl2X10}  +  T  X5 

g2 


An  =  ° 


12 


0 


(4.17) 
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with 


A.(tf)  =  0  i=l ,2 , ,10 


All(tf)  2wllXll(tf) 


X12('tf')  1 


The 


components  of  the  gradient  vector  are  given  by: 


\1  =  -  rl(Bl  X1^3)(kPlBlVkPlx3+kIlV 


+  Y~7  (61x1+X3)A4 


gl 


Hkp2  ‘  r2(B2X2  X3)(kp262X2  kP2X3+kl2X10} 


+  T  .  (e2X2  X3)A5 


g2 


Hk  "  rlX9(kPl6lXl+kP2X3+kIlX9)  +  T  X9AA 

ii  gl 


H, 


12 


r2X10(kP262X2  kP2X3+kl2X10)  +  T  „  x 


g2 


10A5 


(A. 18) 


The  system  canonic  equations  have  been  programmed  for  the  computer 
with  system  parameters  [15] : 


M  =  0.0A  G  =  0.01  T  =  0.5  T  =  0.5  E  =  0.03 
11  gl  tl  1 


M_  =  0.03  G_  =  0 . 008  T  =  1.2  T  =  0.5  D  =  0.5  E0 

2  2  g2  t2  2  2 


0.013 
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and  performance  index  weighting  coefficients: 


i  =  1,2 


8 


r  =  1 
i 


i  =  1,2 


1ft  the  case  that  area  1  (steam)  is  subjected  to  a  different  step 
load  change  and  area  2  (hydro)  is  disturbance-free,  it  has  been  found  that 
the  optimal  control  parameters  of  the  steam  area,  kpl,  k  n,  are  independent 
of  the  magnitude  and  the  sign  of  the  area  disturbance.  For  example,  by 
employing  the  gradient  as  determined  by  equation  1  and  3  of  (4.18)  over 

[tQtf]  =  [0,  20  sec] 

for  a  load  disturbance  of  ±  0.005  p.u.,  conjugate-gradient-descent  has 
been  accomplished  in  5  iterations  (15  trials).  The  initial  guess  of  the 

i  -9 

control  parameter  vector  and  the  step  size  were  taken  equal  to  0  and  10 
respectively.  The  optimal  gain  parameters  have  been  found  to  be  as  follows 


k  =  0.513 

and  the  corresponding  cost 

J  =  0.01436 


with  a  gradient  norm 
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H. 


K 


=  0.2827  x  10 


1 


in  which 


The  augmented  dynamic  system  eigenvalues  are  given  as  follows 
+  0.01485  ±  j  1.82,  -  3.303  ±  j  0.097,  2.8668 

-  0.311  ±  j  0.8418,  and  -  0.641  ±  j  0.068 

Consequently,  the  optimal  P— I  controller  results  in  a  non-asymptot ically 
stable  system.  However,  by  obliterating  the  proportional  part  and  applying 
the  optimization  procedure  with 


k. 


PI 


0 


and 


H 


bi 


0 


the  optimal  integral  control  parameter  obtained  were 


k 


0.5373 


II 


and  the  corresponding  cost 
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J  =  0.01778 

with  a  gradient  norm 

H  | |  =  0.1  x  10"5 
kIl 

The  augmented  dynamic  system  eigenvalues  are  given  as  follows: 

-  0.04996  ±  j  1.751,  -  0.2194  ±  j  7725,  -  3.368, 

-  0.7229  ±  j  0.29145,  -  3.095,  and  -  2.9016 

It  was  found  that  the  optimal  control  parameter  k  ^  is  independent  of 
the  load  change.  Fig.  (6.12)  shows  the  system  transient  and  steady  state 
performance  in  the  case  of  subjecting  area  1  to  a  load  change  of  -  0.005  p.u. 

In  the  case  area  2  (hydro)  is  subjected  to  variable  step  load  changing 
and  area  1  (steam)  regulators  are  inoperative,  it  was  also  found  that  the 

hydro  area  optimal  control  parameters  are  independent  of  the  load  changes. 

If  area  2  is  undergoing,  for  instance,  a  load  disturbance  of  ±  0.005  p.u., 
then  by  using  equation  2  and  4  in  (4.18)  over  a  time  interval 

[tQ,tf]  =  [0,  20  sec.] 

the  conjugate-gradient-descent  was  completed  in  4  iterations  (15  trials) . 

The  initial  guess  of  the  control  parameter  vector  and  the  step  size  were 

_o 

taken  equal  to  0  and  10  respectively.  The  optimal  control  parameters  of 
the  hydro  area  are  as  follows: 


« 


;■ 


TX.  Pt**R  DEV, 
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Fig.  (6-a)  Frequency  deviation-time 


Fig.  (6— b )  Tie-line  power  deviation-time 
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Fig.  (6-c)  Output  power  deviation-time 

System  dynamic  response  of  two  mixed  area  inter¬ 
connected  power  system  obtained  from  non-iteractive 
governor-supplementary  regulator , 


AL1 

=  -0 . 005  p .u . , 

U;L  =  -0.537 

AL2 

=  0.0  p .u .  , 

U2  =  0.0 
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kp2  =  0.6629 

k  =  0.348 

and  the  corresponding  cost 

J  =  0.0206 

with  a  gradient  norm 

=  0.54  x  10~4 

in  which 

K2  *  (kP2’  kI2) 

The  system  augmented  eigenvalues  are  as  follows 

-  3.264  ±  j  0.4596,  -  2.817,  -  1.086, 

-  0.3058,  -  0.2  ±  j.882,  and  -  0.1054  ±  j  1.681 

Fig.  (7)  shows  the  system  transient  and  steady  -state  performance  in  the 

case  of  area  2  being  subjected  to  -  0.005  p.u.  load  changing. 


TJL  POUER  DEV. 
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Fig.  (7-a)  Frequency  deviation-time 


Fig 


(7-b)  Tie-line  power  deviation-time 


« 
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Fig.  (7-c)  Output  power  deviation-time 

System  dynamic  response  of  two  mixed  area  inter¬ 
connected  power  system  obtained  from  non-interactive 
governor-supplementary  regulator , 

AL^  =  0.0  p.u.,  u^=0 

t 

AL^  =  -0.005  p.u.,  u^  =  0.6624  ace^-O.SAS  J  ace9dt 
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The  optimal  interactive  governor-supplementary  regulator 
In  this  mode  of  operation,  the  area’s  supplementary  regulators  operate 
simultaneously  with  the  areas'  governors.  Accordingly,  the  problem  of 
achieving  the  nonintervention  principle  of  the  control  action  cannot  be 
accomplished  by  choosing  the  areas  frequency  bias  setting  equal  to  their 
natural  governing  characteristics  as  in  the  case  of  the  noninteractive 
governor-supplementary  regulator  systems,  since  the  area  control  error  of 
remote  areas  are  not  zero.  Fortunately,  the  sign  of  the  frequency  and 
the  tie-line  power  deviations  may  serve  as  a  criterion  to  show  that  whether 
the  areas  of  an  interconnected  power  system  are  in  need  or  not  to  adjust 
their  own  generation.  If  the  sign  of  the  frequency  and  the  tie-line  power 
deviations  is  the  same  in  one  area,  then  this  area  is  being  subjected  to 
a  local  load  disturbance  and  is  consequently  in  need  of  a  control  action 
to  accommodate  the  load  change.  If  not,  then;  ideally,  the  area  is  not 
locally  disturbed  and  therefore  there  is  no  need  for  its  supplementary 
regulator  to  take  an  action  to  correct  the  mismatch  in  the  prevailing 
load-generation  in  the  remote  areas.  That  is, the  burden  of  the  system 
regulating  responsibilities  are  carried  by  the  disturbed  area;  provided  that 
the  disturbed  area  can  fully  accommodate  its  load  disturbances. 

Basically,  the  problem  of  finding  the  optimal  control  parameters  in  the 
case  of  employing  an  interactive  governor-supplementary  regulator  is  the 
same  as  in  the  case  of  the  noninteractive  one.  The  system  state  and  co¬ 
state  equations  are  the  same  as  in  section  (4.2.1),  except  that  the  initial 
conditions  of  the  state  equations  are  given  as  follows: 


W  = 


0 


i  =  1,2,  ...  ,11 


(4.19) 
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Three  cases  have  been  considered  in  the  problem  of  finding  the 
optimal  control  parameters  of  an  interactive  supplementary  regulator, 
i  A  sudden  loss  or  gain  of  load  at  one  area 

When  one  area, for  example, area  1,  of  the  two  area  interconnected 
power  system  is  being  subjected  to  a  step  load  change  (+  ALp) ,  it  is 
assumed  that  the  disturbed  area  supplementary  regulator  will  take  all 
the  regulating  responsibilities  while  the  remote  area  supplementary  regulator 
is  not  operating;  that  is, 

tf 

U1  =  -  kpl  ace1  -  kn  /  ace1  dt 


in  which 


ace 


1 


Two  schemes  have  been  employed  to  design  the  area  1  supplementary 
regulator : 

(a)  Designing  an  optimal  supplementary  regulator  with  frequency 
bias  setting  equal  to  the  area  natural  governing  characteristic. 

(b)  Designing  an  optimal  supplementary  regulator  with  optimal 
frequency  bias  setting. 

In  the  first  scheme,  the  bias  setting  is  given  by  the  following 
identity 


Bi  ■ 


Ei  +  Gi 


(A. 21) 


, 
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In  this  case  the  control  parameter  vector  is  defined  by 


(4.22) 


The  optimal  control  parameters  were  found  to  be  independent  of  the  load 
changing.  For  example,  in  the  case  that  area  1  is  being  subjected  to  a 
load  change,  ±  0.005  p.u.  the  conjugate-gradient-descent  was  completed  in 


4  iterations  (14  trials) .  The  initial  guess  of  the  control  parameter 
vector  and  the  step  size  were  taken  equal  to  _0  and  10"2  respectively.  The 


optimal  control  parameters  are  as  follows 


=  (0.08335,  0.451) 


(4.23) 


and  the  corresponding  cost 

J  =  0.0242 

with  a  gradient  norm 

I | H  ||  =  0.2  x  10~4 
K1 

In  the  second  scheme,  the  problem  of  designing  an  optimal  supplementary 
regulator  with  frequency  bias-setting  is  considered.  The  problem  was 
reduced  to  a  single  parameter  optimization  problem.  The  control  gain 
parameters  and  k  were  taken  from  identity  (4.23).  In  this  case  the 
gradient  is  given  as  follows 


■ 
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V  =  -  rlkPlXl[kpi(BlXl  +  x3)  +  kUx9' 


+ 


gl 


X2  kFX 


^  -  v9 


(A. 24) 


Adopting  the  technique  which  was  developed  in  Chapter  2,  it  was 
found  that  the  optimal  bias  setting  is  independent  of  the  load  changing. 

In  the  case  of  ±  0.005  p.u.  load  changing,  the  initial  guess  of  the 
frequency  bias-setting  and  the  initial  step  size  were  taken  equal  to  0.04 
(natural  governing  characteristic)  and  0.1  respectively.  After  3  iterations 
(20  trials),  the  optimal  bias  setting  is  given  as 


6  =  0.03595 

and  the  corresponding  cost 


J  =  0.02417 


with  a  gradient  norm 

| |H_  ! I  -  0.3133  x  10'5 
61 

Therefore,  the  value  of  the  optimum  frequency  bias-setting  (0.036) 
is  less  than  the  value  of  the  area  natural  governing  characteristic  (0.04). 
It  was  found  that  the  effect  of  the  area  bias-setting  on  the  value  of  the 
cost  functional,  the  optimality  criterion,  is  negligible.  To  show  that, 
as  well  as  the  effect  of  the  system  parameters  on  the  system  performance, 
a  sensitivity  analysis  has  been  conducted.  This  will  be  presented  later. 
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Similarly,  the  problem  of  finding  the  optimal  control  parameters  in 
the  case  that  area  2  is  being  subjected  to  load  changing  (±  AL2>  and 
the  remote  area  supplementary  regulator  is  not  operating;  that  is, 


in  which 


ace2 


or 


ace2 


APt21  +  62 


has  been  studied. 

In  the  case  of  employing  a  tie-line-bias  control  with  frequency 
bias-setting  equal  to  the  area  natural  governing  characteristic;  that 


is 


5 


e2  =  e2  +  G2 


the  optimal  control  parameters  of  area  2  were  found  to  be 


K2  =  (0.607,  0.32) 


in  which 
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For  example,  if  area  2  is  being  subjected  to  a  load  changing  of  ±  0.005  p.u. 
the  corresponding  cost  is  given  by 

J  =  0.0302 

with  a  gradient  norm 

|  |H  I  |  =  0.38639  x  10"A 
K2 

The  initial  control  parameter  and  the  initial  step  size  were  taken  to  be 
jO  and  10  ^  respectively. 

In  the  case  of  designing  the  area  2  tie-line  bias  controller  with 
optimum  bias-setting,  the  gradient  of  the  Hamiltonian  with  respect  to 
the  bias  setting  is  given  by 


-  x^A 


2  10 


(A. 26) 


Using  the  gradient  as  determined  in  equation  (4.26)  with  the  control  gain 
parameters  given  as 


T 


(0.607,  0.32) 


over 
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the  conjugate  gradient  descent  was  accomplished  with  the  program  described 
in  Chapter  2.  The  initial  guess  of  the  frequency  bias-setting  was  taken 
equal  to  2  natural  governing  characteristic  (0.021).  After  7  iterations 
(10  trials) ,  the  optimal  frequency  bias-setting  was  found  to  be 

B2  =  0.0148 

with  a  corresponding  cost  corresponding  to  dt  0.005  p.u.,  load  disturbance 

J  =  0.029985 


and  a  gradient  norm 


0.141  x  10 


-15 


Once  more,  it  can  be  seen  that  the  effect  of  employing  a  tie-line  bias 
controller  with  optimum  frequency  bias  setting  produces  a  negligible 
effect  in  the  system  optimality  criterion,  J.  For  example,  in  the  case 
of  ±  0.005  p.u.,  load  change,  the  loss  of  optimality,  e, 


E  =  ([J 


-  J 


]  /  J 


opt . £  “non  opt . J  '  “non  opt. 8, 


}  x  100 


is  equal  to  0.712  percent. 

Fig.  (8)  shows  the  system  performance  when  area  1  is  being  subjected 
to  -0.005  p.u.,  load  change  and  the  area  2  supplementary  regulator  is  not 


* 


* 


Ti-  PCKCH  DEV.  ANQ.  FTO.  DEV. 

o*o  o*4  o.m  o.it  r0*0* .  -«*<n  t.oo  o 
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Fig.  (8-a)  Frequency  deviation-time 


Fig.  (8-b)  Tie-line  power  deviation 


0/P  POItR  DRV, 


96 


Fig.  (8-c)  Output  power  deviation-time 

System  dynamic  response  of  two  mixed  area  system 
obtained  from  interactive  governor-supplementary 
regulator , 


AL^  =  -  0 .005  p .u . , 


AL 


5 


u,  =  -0.0834  ace.,-0.451  J  ^ace  dt 


0 


U2  = 


0.0 


0 
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operating.  Fig.  (9)  shows  the  same  when  area  2  is  being  subjected  to  a 

0.005  p.u.,  load  change  and  the  area  1  supplementary  regulator  is  not 
operating . 

ii  Loss  of  a  tie— line  between  two  areas 

This  is  equivalent  to  a  sudden  load  increase  in  one  area  and  a 
sudden  load  decrease  in  the  other.  The  system  parameters  before  and  after 
the  removal  of  the  tie-line  are  the  same  in  the  case  of  an  d.c.  line 
and  not  the  same  in  the  case  that  the  line  is  a.c. 

In  this  case,  both  areas’  supplementary  regulator  should  come  into 
action  simultaneously.  Accordingly,  the  control  law  of  the  IMAIPS  can  be 
written  as  follows 


U1  =  -  kP1  acex  -  kn  /  ace  jdt 

t0 


u2  '  "  kp2  ace2  '  ki2  /  ace2dt 

ko 


in  which 


(A. 27) 


a“l  =  APtl2  +  BlA“l 


ace,  = 


iPt21  + 


f3_  Aw, 


(A. 28) 


For  example,  if  the  system  is  subjected  to  a  load  disturbance  [8], 


AL  =  -  0.005 


AL2  =  -I-  0.005 


[15] 


4 


Fig.  (9-a)  Frequency  deviation-time 


Fig.  (9— b )  Tie-line  power  deviation 


(VP  POUER  DEV. 
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Fig.  (9-c)  Output  power  deviation 

System  dynamic  response  of  two  mixed  area  system 
obtained  from  interactive  governor-supplementary 
regulator , 

AL^  =0.0  ,  u^  =  0 

tf 

AL_  =  -0.005  p.u.,  u9  =  -.706  ace?-0.32  J  ace9dt 
1  A  0 
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which  may  represent  a  sustained  removal  of  one  of  the  tie-lines  connecting 


area  1  to  area  2  . 


The  solution  of  the  canonic  equations  given  by  equations  (4.16), 
and  (4.18)  subject  to  the  boundary  conditions 


i  =  1,2 


12 


yields  the  control  parameter  vector 


(4.29) 


By  using  the  technique  developed  in  Chapter  2,  the  optimum  value  of  the 
control  parameter  vector  was  found  to  be 

K  =  col.  (-.4484,  -.0025,  .3472,  .442) 

and  the  value  of  the  optimum  cost  function 

J  =  0.0192 

with  a  gradient  norm 


H. 


0.971  x  10 


-4 


K 


The  initial  guess  of  the  control  parameter  and  the  initial  step  size  have 


been  taken  equal  to  0  and  0.1  respectively. 

In  the  case  of  employing  an  integral  control  action;  that  is, 


‘ 
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K 


col . 


kI2> 


(4.30) 


the  optimal  control  parameters  were  found  to  be 


K  =  col.  (0.255,  0.27) 


and  the  value  of  the  optimum  cost  function 


J  =  0.0276 


In  the  case  of  applying  the  control  parameters  obtained  in  [15] 
which  are  given  as  follows 


k  =  0.09 


0.4 


The  corresponding  cost 


J  =  0.06946 


which  is  almost  2.5  times  as  the  proposed  controller.  Figure  (10)  shows 
a  comparison  between  the  system  performance  in  both  cases. 


. 
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Fig.  (10-b)  Area  2  frequency  deviation-time 
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Fig.  (10-c)  Tie-line  power  deviation-time 


Fig.  (10-d)  Output  power  deviation-time 


RR^n  oorrmou 
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Fig.  (10-e)  Areas'  control-time 
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Fig.  (10— f )  J  acedt-time 

System  dynamic  response  of  two  mixed  area 

system  obtained  from  the  proposed  controller 

u  =  -  0.0255  I  ace.dt 
1  *  1 

u0  =  -  0.27  J  ace^dt 

and  the  controller  advised  in  [15] 

=  -0.09  /  ace^dt 

=  -0.4  J  ace^dt 

AL^  =  -  0.05  p.u.,  AL^  =  +  0.005  p.u. 
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i i i  A  sudden  loss  of  generation  at  one  area 

This  is  equivalent  to  case  I  except  that  the  system  parameters  before 
and  after  the  disturbance  are  not  the  same.  The  value  of  the  disturbed 
area  inertia  constant  is  lower  after  the  occurrence  of  the  disturbance 
for  instance,  if  area  1  undergoes  loss  of  -  p.u.,  of  its  generation, 

the  new  value  of  Mp  is 

=  Mi  (1  -  AV  (4.31) 

in  which 

n 

Vil  ~  area  1  inertia  constant  after  the  occurrence  of  the  disturbance 

^2.  —  nrea  1  inertia  constant  before  the  occurrence  of  the  disturbance 

It  was  found  that  values  of  the  optimal  control  parameters  depend  upon 
the  values  of  the  areas’  inertia  constant.  Table  (1)  shows  the  relationship 
between  area  1  inertia  constant  and  the  optimal  control  parameters.  Table 
(2)  shows  the  same  for  area  2.  However,  the  optimal  control  parameters 
which  were  obtained  in  the  case  when  considering  that  the  areas’  inertia 
are  constant  results  in  a  stable  system  with  a  slight  loss  of  optimality 
compared  to  the  one  in  which  the  areas’  inertia  were  considered  variable. 

For  example,  in  the  case  that  area  1  is  being  subjected  to  a  0.1  p.u. 
loss  of  its  generation  and  area  1  is  employing  a  supplementary  regulator 
with  gain  parameters 

kpl  =  0.0833 


0.451 


-  ;  :  t  be 


Generation 

loss 

Inertia 

constant 

kpi 

kn 

-0.005 

0.0398 

0.0833 

0.4510 

-0.01 

0.0396 

0.0823 

0.4511 

-0.02 

0.0392 

0.0792 

0.4513 

-0.04 

0.0384 

0.0755 

0.4514 

-0.06 

0.0376 

0.0721 

0.4512 

-0.08 

0.0368 

0.0694 

0.4506 

-0.1 

0.036 

0.0673 

0.450 

Table  (1) 

Steam  area  optimum  control  parameters 


versus  the  area  inertia 


'  '  . 
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Generation 

loss 

Inertia 

constant 

S>2 

kl2 

-0.005 

0.02985 

0.605 

0.321 

-0.01 

0.0297 

0.603 

0.321 

-0.02 

0.0294 

0.500 

0.322 

-0.04 

0.0288 

0.592 

0.324 

-0.06 

0.0282 

0.585 

0.326 

-0.08 

0.0276 

0.57 

0.328 

-0.10 

0.027 

0.57 

0.33 

Table  (2) 

Hydro  area  optimum  control  parameters 
versus  the  area  inertia 


, 
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which  resulted  from  solving  the  system  canonic  equations  with  the 

assumption  that  the  system  parameters  before  and  after  the  disturbance 
the  same;  that  is 

=  M  =  0.04 

this  result  in  the  cost  function 
J1  =  0.099361  x  10+1 

when  it  is  applied  to  the  system  with  =  0.036.  In  the  case  that  the 
system  parameters  before  and  after  are  not  considered  the  same;  that  is 

M  =0.04 
M^1  =  0.036 

The  optimal  control  parameters  are  given  as 

k  =  0.0673 


and  the  value  of  the  optimum  cost  function 
J9  =  0.9928  x  10+1 

4— 


are 


Hence,  in  the  case  of  employing  a  control  strategy  which  neglects  the 
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changes  in  the  systems’ 


inertia  results  in  a  loss  of  optimality 


£ 


VJ2 


X  100  =  0.0815 


percent 


Therefore,  the  assumption  that  the  systems’  inertia  are  equal  before  and 
after  a  generation  disturbance  is  satisfactory. 


4.3  Sensitivity  analysis  of  the  LFC  problem 

Sensitivity  considerations  of  dynamical  systems  are  very  important 
m  designing  control  systems.  In  designing  a  control  system  one  must  first 
develop  a  mathematical  model  which  represents  the  physical  behavior  of 
these  systems.  There  is  always  a  discrepancy  between  the  physical 
reality  and  the  mathematical  model.  This  is  due  to;  inexact  identification 

of  the  systems  parameters,  unpredictable  changes  in  the  systems  environment, 
etc . 

Sensitivity  analysis  is  of  a  particular  importance  in  the  case  of 
employing  modern  control  theory  to  design  systems  with  prescribed  behavior. 
The  results  obtained  employing  modern  control  theory  are  useless  in  practice, 
if  the  mathematical  model  used, prove  to  be  very  sensitive  to  parameter 
variations.  Therefore,  it  is  very  important  to  perform  sensitivity  analysis 
whenever  modern  control  theory  is  adopted. 

A  sensitivity  analysis  has  been  performed  on  the  problem  of  LFC  of 
two  mixed  area  interconnected  power  system  to  study  the  system  performance 
to  the  system  parameters  as  well  as  to  the  control  parameter  variations. 

The  sensitivity  analysis  of  a  dynamic  system  can  be  defined  as  follows: 

Given  a  dynamic  system  described  by  the  following  ordinary  differential 


equat ion 


‘ 


Ill 


Xj_  =  A1(K)X1+riAL,X1 


0 


(A. 32) 


in  which 


X^  is  an  nxl  state  vector 

K  is  an  mxl  parameter  vector 

the  problem  posed  is  to  find  the  sensitivity  of  the  system  to  the  parameter 
variations,  AK.  One  possible  way  to  find  the  system  sensitivity  to  the 
parameter  variation  AK  is  to  find  the  sensitivity  of  a  cost  functional 

1  ^  T 

J  =  2  J  XlQlXldt  (4.33) 

t 

0 

with  respect  to  a  variation  in  the  parameter  vector  AK. 

To  handle  this  problem,  one  can  introduce  a  new  variable 

*0  =  X1Q1X1  x0(t0)=  0  (4*34) 

and  by  employing  the  trajectory  sensitivity  analysis  of  dynamic  systems 
one  can  find  the  sensitivity  of  x  (tf)  with  respect  to  the  parameter  K. 

Then,  define  a  new  variable  X  as  follows 

X  =  col.  (X1,  xq)  (4.35) 


The  augmented  dynamic  system  can  be  written  as 
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X  =  A(K)  X  +  TL,  x(t0)  =  0 


(4.36) 


By  taking  the  partial  derivative  of  equation  (4.36)  with  respect  to  one 
obtains 


ax  _  9A(K)  3X  9A(K) X  9X(t0)  _ 

ak  ax  ak.  ak.  »  Ik 

1  1  1  -i 


or 


=  9A(K)  ^  3A(K) X  .  .  „ 

i  ax  Si  +  3k.  ’  Si^0  “  °»  1  ■  1,2, ...,m  (4.37) 


in  which 


si(t) 


is  an  [ (n+1) xl ]  sensitivity  vector, 


Therefore,  the  cost  functional  sensitivity  with  respect  to  the  elements 
of  K  is  given  as 


■aiT  =  Si(tf)  i=l,2,...,ra  (4.38) 

i 

Equation  (4.38)  can  be  easily  calculated  by  the  forward  integration 
of  equation  (4.36)  and  equation  (4.37). 

In  the  case  that  area  1  is  being  subjected  to  -0.005  p.u.  load  change, 
and  the  nonintervention  principle  of  the  areas  supplementary  regulator  is 
satisfied,  the  cost  functional  sensitivity  with  respect  to  area  1  control 
parameter  variations  is  as  follows 
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S  =  -  .235  x  10 
kPl 

S,  =  -  .426  x  10  2 
II 

S  =  -  0.3883  x  10"3 
B1 

For  a  0.05  change  in  the  control  parameters  kpl>  Is.  and  B  the  percentage 
change  in  the  value  of  the  cost  function  is  -  0.004,  -  0.397,  and  -  0.003 
respectively.  This  will  emphasize  the  effect  of  the  integral  control  parameter 
on  the  system  performance. 

The  sensitivity  of  the  cost  function  with  respect  to  the  area  1 
parameter  variations  are 

S  =  -  .469  x  10-2 
gl 

S  =  +  .  145  x  10"3 
tl 

=  -  .  415  x  10"2 

E1 

For  example,  5  percent  change  in  T  ^ ,  T^,  and  results  in  -  0.484, 

+  0.015,  and  -  0.026  percentage  change  in  the  value  of  the  cost  function. 

This  shows  the  prominent  effect  of  the  governor  time  constant  on  the 
system  performance  compared  to  that  of  the  turbine  time  constant  and  the 
speed  regulation.  Therefore,  more  consideration  should  be  given  to  the 
process  of  identifying  the  areas'  governor  time  constant  than  to  that  of 
the  other  two  parameters. 
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Similar  results  were  obtained  from  area  2  sensitivity  analysis. 

For  example,  when  area  2  is  being  subjected  to  -  0.005  p.u.  load  change, 

the  cost  functional  sensitivity  with  respect  to  area  control  parameters 
were  found  to  be 

S,  =  -  0.723  X  10~2 

h>2 

S.  =  -  0.10  x  lcf1 
kI2 

SD  -  -  0.733  x  10~2 
82 

e.g.  for  a  5  percent  change  in  the  area  control  parameters  ,  k^  >  and 
^2 »  there  is  a  0.0725,  0.532,  and  0.0255  percentage  change  in  the  value  of 
the  cost  functional. 

The  sensitivity  of  the  cost  functional  with  respect  to  area  2 
parameter  variations  are  as  follows 

S  =  -  0.281  x  10~2 
g2 


S  =  -  0.25  x  10 
t2 

S  =  -  0.421  x  10-3 
2 

S  =  -  0.12115  x  10 
E2 


which  again  show  that  the  area  2  governor  time  constant  has  a  conspicuous 
influence  on  the  system  performance. 
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4.4  Nonlinear  load  frequency  control  of  interconnected  power  system 

The  problem  of  load  frequency  control  (LFC)  exhibits  many  nonlinearities. 
The  most  significant  of  these  arise  from  the  tie-lines  connecting  the  areas 
of  an  interconnected  system,  and  the  dead-band  of  the  governors.  The 
former  is  overriding  when  the  system  is  subjected  to  large  disturbances 
which  produce  large  power  angle  deviations  while  the  latter  is  prominent 
for  small  disturbances.  In  this  section,  the  tie-line  nonlinearity  is 
considered . 

The  tie-line  nonlinearity  arises  from  the  relation  between  the  tie-line 
power  and  the  power  angle  deviation.  The  tie— line  power  deviation  of  area 
i  of  an  n-area  interconnected  power  system  is  given  by 

n 

Ap  •  =  I  T  sin  6  (cos A <5 .  .  -1)  +  T  .  .  cos6° .  sinA6  .  . 
tx  j=l  1J  13  ij  1 3 

(4.39) 

Miniesy  and  Bohn  [16]  have  proposed  a  two  level  controller  for  the 
problem  of  nonlinear  load  frequency  control.  The  first  level  is  a  local 
feedback  control  for  each  area  of  the  interconnected  system.  This  control 
is  a  function  of  its  own  area  state  variables.  The  second  level  is  an 
intervention  open  loop  which  is  utilized  to  compensate  for  neglecting  the 
coupling  state  variable  between  the  areas  of  the  interconnected  system, 
and  the  nonlinearities  due  to  the  tie-lines.  The  idea  of  employing  a 
multi-level  control  strategy  in  the  area  of  load  frequency  control  is  not 
adequate.  First,  it  increases  the  controller  design  complexity.  Second; 
it  is  not  used  in  practice.  Doraiswami  [17]  has  considered  the  problem  of 
LFC  with  a  tie-line  nonlinearity  from  the  stochastic  viewpoint.  An  observer 
was  employed  to  implement  the  control  law.  Once  an  observer  is  introduced 
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into  the  system,  the  cost  increases  and  the  control  is  no  longer 
optimal.  Therefore,  it  is  necessary  to  design  a  one  level  optimal 
controller , and  the  control  law  must  only  be  a  function  of  the  measurable 
states . 


A. 4.1  The  mathematical  model 

The  dynamic  nonlinear  model  of  steam-hydro  interconnected  power 
systen  (SHIPS),  is  the  same  as  was  considered  in  the  preceeding  sections, 
except  that  the  state  variable  x^ ;  which  was  representing  the  tie-line 
power  deviation,  is  representing  the  power  angle  deviation.  Therefore, 
the  dynamic  model  of  SHIPS,  taking  into  account  the  tie-line  nonlinearity, 
is  given  by  the  following  system  of  equations 


x 


1 


x 


2 


X1  '  (Acosx3+  Bsinx3-  A)  +  i-  xfe 

G2  1  1 

x2  +  g-  (Acosx3+  Bsinx3-  A)  +  —  xg 


AL 

AL 

ML 


1 


2 


x3  =  xi  "  x2 


X4  = 


1  1  ,1 

T  X1  T  X4  T  U1 


gl 


gl  gl 


2  1  ,1 

x  =  -  — - — —  xc  +  — —  u_ 

5  T  T  5  T  2 

g2  g2  g2 


X ,  =  - —  x.  -  ■= X, 

6  Ttl  4  Ttl  6 


.  _  1  1 

ij.  xc  ~  t  X7 

1 1 2  t2 


.  =  2_  2 
X8  Tt2  X5  (Tt2+  D2)  "7  D2  "8 


x„  -  r—  x, 


(4.40) 
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in  which 


A 


B 


sin 


cos 


6 


6 


0 

12 


0 

12 


The  control  law  of  each  area  is  specified  to  be  proportional— 
plus— integral  the  area  control  error.  In  terms  of  the  systems’ 
state  variable,  up  and  can  be  written  as 

t{ 

U1 - kpl(B1x1+ Acosx3+  Bsinx3-A)  -  k  J  (B.x  +  Acosx 

lC 


+  Bsinx^-  A)dt 


u2  ~  ~  k^2($2x2~  Acosx3~  Bsinx^^ 


A)  J 


0 


(82x2~  Acosx^ 


-  (3sinx^+  A)dt 


(A. 41) 


The  problem  posed  is  to  find  the  optimal  controls  up  and  u2  or 
equivalently  the  optimal  parameters  kpi ,  kp2 ,  k  and  kJ2  which  minimize 
the  cost  functional 


-k  I 


0 


8 

(  l 

i=l 

i^3 


qiXi  +  (Acosx3+  Bsinx^-  A) 


2 

+ 1 

i=l 


r .u . )dt 

l  l 


(4.42) 
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subject  to  the  constraints  given  by  equation (4. 41)  and  equation  (4. 42). 
Also,  the  steam  rate  of  generation  should  be  less  or  equal  Y  p.u. 
mw/ sec . 

To  cast  the  system  dynamic  equations  into  the  familiar  state  space 
form,  one  might  introduce  new  variables  x^  and  x^q  such  that 

Xg  =  6^^  +  Acosx^  +  Bsinx^  -  A  (4.43) 


and 


*10  =  ^2*2  ”  AcosX3  “  Bsinx^  +  A 

to  formulate  the  control  law  into  the  proportional  form.  Consequently, 
u^  and  u^  can  be  written  as  follows 


ul  =  ~  kpi^ixi+  Acosx3+  Bsinx^-  A)  -  ki;1xg 


u2  =  “  kp2^2X2~  Acosx3~  Bsinx3+  A)  -  kI2  x  Q  (4.44) 


The  inequality  constraints  due  to  the  steam  area  rate  of  generation 
limit  can  be  expressed  mathematically  as  follows 


(4.45) 


To  convert  equation  (4.45)  into  an  equality  constraint,  another  new 
variable  x^  is  introduced 


4 
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-  *>  Hl(? 


tl 


(A. 46) 


where 


H^ci) 


if  a  <  0 


if  a  >  0 


The  problem  of  finding  the  optimal  control  parameters  (k  ,  k  , 

K  PI’  P2* 

^11’  anc^  ^12 ^  system;  given  by  equation  (4.41),  equation 

(4.42),  equation  (4.44)  and  equation  (4.47),  using  the  technique  which 
has  been  introduced  in  Chapter  2  failed  to  converge.  Therefore, 
another  approach  was  adopted  by  reformulating  the  system  dynamics 
into  a  quadratic  form.  This  was  done  by  introducing  a  nonlinear 
transformation  [l8],  on  the  angular  frequency  deviation  (x^).  This 
nonlinear  transformation  is  defined  by 


x12  =  sinx3 


and , 


then, 


define  a  pseduo-variable, 


=  /l  - 


x 


2 

12 


(4.47) 


(4.48) 


For  convenient  notation,  a  new  variable  y  is  introduced  such  that 

y.  -  x.  i  =  1,2 

J  i  i 

y.  =  xj+1  j  =  3, 4, 5,..., 10 


-• 
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then,  the  system  dynamic  equations  can  be  written  as  follows 


>1  ■ 

G1  1 

=  -  yl  -  <A  yi2  +  B  yn  -  A>  +  ^  y5  -  ^  AL1 

y2  = 

G2  1 

■  -  s;  y2  +  (A  y12  +  B  yn  -  A>  +  RJ  y7  -  al2 

E1  1  ,1 

T  yi  T  y3  T  U1 

gl  gl  gl  1 

• 

y4  = 

E2  1  1 

T  y2  T  y4  +  T  u? 

2  g2  g2  Z 

• 

y5  = 

1  1 

Ttl  ^  Ttl  "5 

• 

y6  = 

Tt2  y4  -  Tt2  y6 

II 

2  ...  2  2 

Tt2  ^  (Tt2+  V  Yfi  D2  Y? 

^8 

6lyl  +  yl2  +  B  yll  "  A) 

y9 

e2y2  ~  (A  y12  +  B  yn  "  A) 

yio 

,1  1  ,  n2  „  ,1  1 

-  (  y0  ~  T  Yc  V)  H-1  (T  y-5  ”  T  ys  " 

tl  J  tl  tl  J  tl 

yll 

"  y12(yl  '  y2> 

yl2 

=  yn(y2  -  y2) 

(A. 49) 
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and 


U1  kPl(6lyl  +  A  y12  +  B  ^11  -  A>  -  knyg  (4.50) 

U2  =  •  kP2(B2y2  *  A  y12  *  B  yll  +  A>  -  kI2y9 

To  handle  a  criterion  of  the  form  (4.43),  one  can  augment  the 
system  differential  equations  by  an  additional  equation  as  follows 

1  9  2 

y13  =  2  { Vi  +  VA  yi2  +  B  yn  - 

+[kPl  (Blyl  +  A  y12  +  B  Vn  *  A)+  knyg]2 

+[kP2  (B2y2  "  A  yl2  -  B  yll  +  A)2+  ki2y9]2}  <4.51) 

then,  by  defining  a  Hamiltonian 

13 

H  =  l  X.y. 
i=l 

and  applying  Pontryagin’s  minimum  principle,  one  can  get  the  system  co¬ 
state  equations 


Ai  qiyi  "  ri6ikpi^kpi(eiyi  +  A  y12  +  B  vn  ‘  A)  +  kny8> 


G  E  6 

+  mT  xi  +  (fT  +  T~  kPl)  X3  "  eiX8  "  yi2  A11  +  yn  A 12 
1  gl  gl 


•  . 
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A2  "  q2y2  "  r262kP2{kP2(62y2  '  A  yl2  "  B  yll  +  A)  +  kT2y9* 


12  "  '11 


'12^9 


G2  E2  B2 

t>2*2+  <Tg2  +  Tg2  kp2)  *4  '  B2X9  +  y12XH  ‘  yllX12 


A.  =  - 


q3y3  +  T 


A3  "  T 


,  _  2_  1_ 

A5  T.  ,  kT.  ,  y3 


gl  tl  xtl  xtl  J  xtl 


T.  ,  y5  ^H1A10 


A  =  - 


q4y4  +  7 


A,  - 


T  „  A6  +  T_  A7 


g2  t2  t2 


A5  =  *  V5  ”  M2  A1  +  AS  +  y3  -  T  y5  "  ^')H1A10 


A6  =  -  V6^Vt^A7 


A7  =  -  q7y7  -  k~2  A2  +t2  A7 


A8  =  -  Vs  -  rlkIl{kPl(Vl  +  A  y12  +  B  y12  -  A)  +  kny  } 


+  FT  kn  A3 

gl 


x9  =  “  99Y9  -  r2kI2{kp2(62y2  "  A  Yl2  "  B  y12  +  A>  +  ki2y9} 


+  tT  kI2  A4 

g2 


‘10 


0 
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X11  =  -  V(A  yl2  +  B  yu  *  A)  -  rpkppA{kpp(Spyp  +  A  y 


+  B  y  u  -  A)  +  kn  yg} 


r2kP2A^kP2  ^2y2  A  y 


12 


B  yn  +  A)  +  kI2y9> 


+  —  i  .  A_  . 
M1  A1  M2  X2 


+  t  kpiA3  t  kP2A4  -  A  A8  +  A  *9  -  (y2  -  yx)A 
gl  g2 


12 


A12  «  -  qTB(A  yu  +  B  y^  -  A)  -  r!kPlB{kPl  (61^1  +  A  yl2 


+  By  —  A )  +  k  y  } 
>11  '  Ily8 


+  r2kp2Blkp2(62y2 


12 


■Byn+A)  +  kl2y9} 


M1  A1  M2  A2 


+ 


B 


gl 


kPlX3  * 


B 


g2 


kP2\ 


-  B 


X8  +  B 


x9 '  (yry2)xn 


A 


13 


=  0.0 


(4.52) 
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and  the  gradient  components  are  given  by 


Hkpi  ‘  rl(6lyl  +  A  y12  +  B  yn  “  A)  +  A  y12  +  B  yu  -  A) 


+  kIl  y8}  +  (6lyi  +  A  y12  +  B  yu  -  A) A. 


~  r2(62y2  “  A  yi2  "  B  yn  +  A)  {(B2y2  "  A  yl2  "  B  yll  +  A) 


P2 


+  kI2y8}  + 


T —  (B0y 


g2 


2J  2 


-  A 


12 


-  B  yn  +  A>X4 


rly8{kPl(Vl  +  A 


12  +  B  yu  -  A)  +  k  y  }  +  —  y  X 

gl 


r2y9(kP2(62y2  A  y12  ~  B  yll  +  A)  +  ki2  y9}  +  T  y9  X4 

g2 


(4.53) 


The  non-intervention  principle  of  che  control  actions  was  the 
strategy  adopted  in  finding  the  optimal  control  parameters.  It  is 
assumed  that  both  speed  regulators  and  governors  are  operating  simultaneous¬ 
ly  in  response  to  the  load  changing  in  its  own  area. 

First  area  1,  the  steam  area,  was  subjected  to  a  different  step 
load  changing.  It  was  assumed  that  the  initial  load  angle  difference 
between  area  1  and  area  2  was  given  as 
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6 j  £  ~  1*0  radian 

As  was  expected,  when  the  system  was  subjected  to  small  load 
disturbances,  the  optimal  control  parameters  were  found  to  be  the  same 
as  in  the  case  of  the  linear  load  frequency  control  problem.  Table  1 
shows  the  optimal  control  parameters;  kpi  and  k  for  different  step 
load  changing. 

It  can  be  seen  from  Table  3  that  the  optimal  control  parameters 
are  not  independent  any  more  of  the  load  changing  and  an  adaptive 
controller  is  needed  to  get  optimal  performance  for  each  load  disturbance. 

Second,  area  2  the  hydro  area,  was  subjected  to  different  step  load 
changing.  It  was  assumed  that  the  initial  operating  condition 
characterized  by  the  load  angle  difference  between  area  1  and  area  2 
was  given  as  follows 

0 

^21  ~  radian 

Similar  results  and  conclusions  were  obtained  as  in  the  case  of  disturbing 
the  steam  area.  These  results  are  summarized  in  Table  4. 

4.4.2  The  effect  of  the  initial  power  angle  on  the 

optimal  control  parameters 

The  conclusion  that  the  areas’  optimal  control  parameters  depend 
on  its  area  load  changing  was  the  motivation  to  study  the  effect  of  the 
initial  power  angle  on  the  areas’  control  parameters. 
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6°2  rad 

0.4 

0.6 

0.8 

1.0 

kpi 

0.13 

0.114 

0.1 

0.1 

kn 

0.636 

0.600 

0.563 

0.499 

Table  5 

Steam  area  optimum  control  parameters 
versus  area  initial  power  angle 


6°^  rad 

0.3 

0.5 

0.7 

0.9 

kP2 

1.109 

1.1025 

1.059 

0.976 

kI2 

0.3193 

0.321 

0.32  7 

0.337 

Table  6 

Hydro  area  optimum  control  parameters 
versus  area  initial  power  angle 


• 

129 


For  example,  if  area  1  is  being  subjected  to  +  0.05  p.u.  load 
disturbance,  the  relation  between  the  initial  load  angle  and  the 
optimal  control  parameters  is  shown  in  Table  5.  Similar  results (shown 

in  Table  6  are  obtained  when  area  2  is  being  subjected  to  a  0.05  p.u. 
load  disturbance. 

It  was  noticed  that  in  the  case  of  neglecting  the  tie-line  non¬ 
linearity  in  the  system  dynamics  that  the  system  is  always  stable  and 
can  accommodate  any  amount  of  load  disturbances.  In  contrast;  in  the 
case  of  taking  the  nonlinearity  of  the  tie-line  into  consideration,  there 
is  a  load  increase  limit  for  each  power  angle.  For  example,  if  6^ 
is  equal  to  1.0  rad,  the  maximum  load  which  area  2  can  fully  accommodate 
is  0.025  p.u.  In  the  case  of  6^  is  equal  to  -1 .  radian,  the  maximum 
load  increase  vrhich  area  2  can  withstand  without  losing  its  stability 
is  0.15  p.u . 

In  conclusion,  the  nonlinear  LFC  control  parameters  depend  on  the 
initial  power  angle  and  the  megnitude  of  the  load  disturbance.  However, 
the  application  of  the  optimal  control  parameters  which  were  obtained  from 
the  linear  analysis  of  the  LFC  to  the  nonlinear  version  of  the  problem 
results  in  an  insignificant  loss  in  the  optimality  performance  index.  For 
example,  if  area  1  is  being  subjected  to  a  0.06  p.u.  load  change  and 
^12  iS  Gqual  to  1,0»  the  value  of  the  cost  function  in  the  case  of 
employing  the  linear  optimal  control  parameters  to  control  the  nonlinear 
system  is 


J  =  3.094 


- 

•  : 


. 
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On  the  other  hand,  in  the  case  of  employing  the  nonlinear  optimal 
control  parameter,  the  value  of  the  cost  function  is 

J2  =  3.053 

The  number 

J1  ”  J2 

e  =  - 7 1  x  100 

J1 

or 

e  =  1.325  percent 

can  serve  as  a  suboptimality  index  for  the  nonlinear  system  under  employing 
the  linear  optimal  control  parameters  (suboptimal  control  parameters)  to  control 
the  nonlinear  system.  Therefore,  e  represents  the  price  we  have  to  pay,  in 
terms  of  the  performance  index,  to  achieve  better  dynamic  performance  of  the 
nonlinear  system.  Hov’ever ,  the  loss  in  the  performance,  in  the  case  of 
employing  the  suboptimal  control  parameters  instead  of  the  optimal  one  is 
insignificant  compared  to  the  complications  that  arise  in  the  design  of  the 
control  system  due  to  the  necessity  of  the  introduction  of  an  observer  in 
the  system  to  identify  the  load  change  and  consequently  the  corresponding 
optimal  control  parameters.  Figure  11  shows  a  comparison  between  the  system 
performance  in  the  case  of  employing  the  optimal  nonlinear  control  parameters 
and  the  linear  one  to  control  the  nonlinear  system. 


. 
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Fig.  (11-a)  Area  1  frequency  deviation-time 


Fig.  (11-b)  Area  2  frequency  deviation -time 
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Fig.  (11-d)  Output  power  deviat ion-time 


O.W  _  -0.0«  -0.04  -O.Ot  0.00  O.Ot 
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Fig.  (11-e)  Area  1  control-time 
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TU"E  IN  SECOTCS 

Fig.  (11-f)  I  ace^t  -  time 


Comparison  between  system  dynamic  response  of  two 
area  mixed  system  obtained  from  the  nonlinear 
optimal  controller  parameters;  optimal  solution, 
and  the  linear  optimal  control  parameters;  the 
suboptimal  solution. 


Chapter  5 

Load  Frequency  Control  With  Governor-Dead  Band 


5.1  Introduction 

The  governor-dead  band  is  the  hysteresis  effect  which  arises  from 
mechanical  friction  and/or  looseness  in  the  moving  parts  of  the  speed- 
governing  system.  That  is,  for  any  given  controlled  valve  position,  there 
may  be  a  change  in  the  control  signal  direction  before  the  valve's  position 
changes,  kirchmayer  [6]  has  conducted  an  experimental  analysis  using 
an  analog  computer. to  study  the  effect  of  the  governor-dead  band  on  the 
system  behavior  when  the  system  is  subjected  to  step  load  changes.  It 
was  shown  that  the  dead-band  effect  is  significant  when  the  system  is 
subjected  to  small  load  changes. 

A  governor  hysteresis  input-output  characteristic  may  be  described 
by  the  following  equation 


x .  -A 

/  ln 

X  =<  X.  +  A 

out  m 

V 

constant 


x . 
m 


x . 
m 


x . 
m 


>  0 


<  0 


-  x 

out 


(5.1) 


<  A 


in  which 


2A  =  dead  band 
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It  is  clear  that  hysteresis  makes  the  system  past  history  very 
important  in  influencing  the  system  future  response.  This  is  in  contrast 
to  the  systems  considered  in  the  preceeding  two  chapters,  in  which  a 
previous  state  has  no  effect  once  the  system  trajectory  has  carried  the 
system  away  from  that  state. 

The  problem  of  finding  the  optimal  controls  for  a  system  which 
incorporates  a  hysteresis  nonlinearity  will  be  handled  in  the  same  manner 
in  which  systems  without  a  hysteresis  nonlinearity  were  handled  in  the 
preceeding  chapters.  The  chief  difference  is  that  the  adjoint  system 
differs  slightly  from  the  one  in  which  the  hysteresis  nonlinearity  is  not 
included.  To  show  how  to  deal  with  systems  having  a  hysteresis  loop  in 
an  optimal  sense,  the  problem  of  LFC  of  a  single  steam  area  will  be 
considered . 

5.2  Load  frequency  control  of  single  steam  area  power  svstem  with 
dead  band 

A  typical  linearized  model  of  a  steam  area  power  svstem  with  a 
governor  dead-band  is  shown  in  Fig.  (12)  [14].  The  control  lav;,  u,  is 
specified  to  be  proportional-plus-integral  of  the  angular  frequency  deviation. 
In  state  space  the  system  dynamics  can  be  written  as  follows 
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the  control  law  can  be  written  as  follows 


u  =  - 


k  xi 

P  1 


kI  J 
0 


x^dt 


(5.3) 


in  which 


x2  -  A 


x2  <  0 


x20=<x.  +  A 


constant 


x  2  >  0 


|x2o-x2[  1  A 


To  cast  the  system  dynamics  into  the  familiar  state  space  form, 
a  new  variable;  x^,  is  introduced 


and  therefore  the  control  law  can  be  written  as 


u  =  -  k  x,  -  k  x ^ 
pi  15 


(5.5) 


k  and 
P 


The  problem  posed  is  to  find  the  optimal 
k^.,  which  minimize  the  cost  functional 


control 


parameters 


1  ib  ?  7 

=  J  J  {  I  Vi  +  r(kDxi+  kix5)  }dt 
0  i=l  p 


(5.6) 


subject  to  the  dynamic  constraints  given  by  equation  (5.2).  Here  the 
rate  generation  constraint  is  not  considered  since  the  problem  of  LFC  is 
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studied  for  small  load  disturbances. 

To  handle  a  criterion  of  the  form  (5.6)  a  new  state  variable  is 

def ined 


x 


6 


Vi  +  r(kpV  kix5)2) 


(5.7) 


with 


x6(t0)=  0 

therefore,  the  cost  functional  can  be  written  as  follows 


J  =  W 


(5.8) 


Now,  by  defining  a  Hamiltonian 


H 
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I  X 

i=l 


.x . 

i  i 


(5.9) 


then,  by  applying  Pontryagin’s  minimum  principle,  one  can  get  the  system 
equations.  The  system  state  equations  are  given  by  equations  (5.2),  (5.4) 
and  (5.7).  The  system  co-state  equations  are  found  to  be 
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=  0 


with 


Ai(tf)  =0  i  =  1,2, ... ,6 


in  w7hich 


^  1 


dead  band  is  active 


otherwise 


the  gradient  vector  components  are  given  by 


(5.10b) 
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xlX2 


(5.11) 


Since  the  integrand  in  the  cost  functional  equation,  equation 
(5.6),  and  the  system  dynamics  are  not  explicit  functions  of  t,  we  have 


H  =  0 
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or 


H  -  constant  on  the  optimal  trajectory 


Therefore,  in  the  dead-hand  zone  the  following  condition 


3H  |  3H  | 

Sx  -  8x  '  + 

d  Cd 

should  be  satisfied,  that  is, 


A  .  (t  )  =  A . (t  )  i=l,2, ... ,6 

id  id 


(5.12) 


and 


Hk  (td} 
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Hk  <0 
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in  which 
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t  ,  =  lim  t  ,  +  e 
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x2(t')  =  x2(tj)  =  o 


where 
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that  is,  the  time  at  which  the  dead  band  starts  or  ceases  to  be  active. 
Equations  (5.12)  are  similar  to  the  Weierstrass-Erdman  corner  conditions. 

A  corner  point  of  x(t),  is  a  point  at  which  the  derivative  is  not  uniquely 
defined,  i.e.,  a  point  at  which  dx/dt  poses  a  jump  discontinuity. 

5.3  The  numerical  solution  of  TPBVP  of  systems  incorporating  a  back  lash 

element 

The  problem  of  finding  the  optimal  control  parameters  of  a  P-I 
controller  of  a  dynamic  system  which  includes  a  back  lash  element  has  been 
reduced  to  the  problem  of  solving  a  nonlinear  two— point  boundary  value 
problem  with  multiple  corner  points. 

To  solve  a  two-point  boundary  value  problem  with  multiple  corner 
points,  the  technique  which  has  been  introduced  in  Chapter  2  with  a  slight 
modification  will  be  adopted  This  modification  will  take  place  in  the 
forward  path  of  integration  to  facilitate  the  detection  of  the  corner 
points.  As  the  integrating  process  in  the  forward  direction  proceeds,  two 
values  of  each  x^Ct^.)  and  ^(t^)  j  =  1-2,  i-1,  are  stored.  When  a  new 
point  x^Ct^)  is  being  integrated,  two  problems  have  to  be  solved  before 
proceeding  with  the  integration  process.  These  are: 

1.  detection  of  the  mode  changing  which  may  be  classified 
as  follows: 

i  -  rising  mode  of  operation 
ii  -  falling  mode  of  operation 
iii  -  dead  band  mode  of  operation 


2. 


locating  the  exact  moment  of  mode  transition,  i.e.,  the 


■ 
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moments  of  entering  and  leaving  the  dead  band 
(corner  points). 

Detection  of  the  change  of  the  mode  of  operation  can  be  easily 
done  as  follows 

(i)  in  the  case  of  entering  the  dead  band  mode,  the  sign 

of  the  prevailing  is  opposite  to  the  previous  one; 
that  is, 

sign{x2(ti)  *  (t  )}  =  -  ve 

(ii)  in  the  case  of  leaving  the  dead  band,  the  following 
inequality  should  be  satisfied: 

|x  (t  )  -  x  (t . )  |  >  A 
2  r  2  i  1  — 

in  which 

x2(t  )  =  the  last  reversing  point 

(iii)  in  the  case  of  a  rising  mode  of  operation: 

sign{x2  (t  ) }  =  +  ve 

(iv)  in  the  case  of  a  falling  mode  of  operation 

signlx^ (t^) }  =  -  ve 
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The  problem  of  finding  the  exact  time  of  entering  the  dead  band, 
reversing  can  be  done  by  employing  one  of  the  interpolation 


techniques .  If  at  the  present  point,  t.,  of  integration 

sign{x  (t . )  x0(t.  )  =  -  ve 

2  1  2  l-l 

is  satisfied,  this  means  that  the  system  has  already  passed  the  reversing 

point  (tj.)*  To  determine  t^_ ,  one  can  employ  a  Hermite  third  order 

interpolation  polynomial  to  approximate  x.(t)  through  x^(t  n)  and  x^(t  ). 

2  2  i-2  2  i-1 

The  third  order  Hermite  polynomial  is  given  by  the  following  equation: 


(5.13) 


in  which 
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This  reduces  to 


x2  (t) 


t3  + 


b  t  +  ct  +  d 


(5.14) 


in  which 


^i-rW 


b  "  {3(ti-l+  ti-2)^x2^t±-l^~  x2(ti-2)  + 


(ti-l'  ti-2>[(ti-l+  2ti-2)x2(ti-l)+2ti-l+ 
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d  {(3ti-r  ti-2)ti-2  X2(ti-l)+(ti-i“3ti_2),:i-l  x2(ti-2) 


^i-l”  ti-2^ti-2X2  ^i-l ti-lVti  2^ 


The  reversing  point  (t  )  occurs  when 


dx2(t) 

T7  =  3a  t  +  2bt  +  c  =  0 
dt  r  r 


or 


t  = 
r 


(5.15) 


Knowing  the  reversing  time,  the  program  re-integrates 


If 


x9  (t  ) 
2  r 


x  (t  )  - 
2  r 


<  E 


where  e  is  a  predetermined  limit 

then  x2^tr^  or  x2^tr^  is  accePtec3  as  a  corner  point.  If  not,  x9  will 

be  replaced  by  x?  (t^)  and  the  process  is  repeated  until  the  reversing 

point  is  located  accurately.  Then  x.(t  )  is  stored. 

2  r 

The  problem  of  finding  the  exact  time  of  leaving  or  entering  the 
dead  band  can  be  solved  by  employing  an  inverse  interpolation  process. 

Once  again  the  Hermite  third  order  polynomial  will  be  used  iteratively  to 
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find  the  exact  point  of  entering  or  leaving  the  dead  band.  This  can  be 
done  by  solving  the  following  equation 

x2(tr)  1  2A  =  a  t^  +  b  t*  +  c  td  +  d  (5.16) 

where  td  is  the  time  at  which  the  dead  band  starts  to  be  inactive.  Equation 

(5.16)  can  be  solved  analytically,  but  it  has  been  found  that  the  analytical 

solution  accuracy  is  dependent  on  the  coefficients  a,  b  and  c.  The 

Newton  Raphson  method  has  been  employed  to  solve  equation  (5.16)  instead 

of  the  analytical  method.  As  in  the  case  of  finding  reversing  points,  the 

process  of  finding  t  will  be  iteratively  repeated  until  t  ,  is  determined 

u  d 

accurately . 

Having  solved  the  problems  of  finding  the  reversing,  entering,  and 
leaving  points,  one  can  solve  the  two— point  boundary  value  problem  given 
by  equation  (5.2)  and  (5.10)  using  the  algorithm  of  chapter  2  to  get  the 
optimal  control  parameters. 

5.4  The  optimal  solution 

The  system  was  subjected  to  different  load  disturbances.  It  was 
found  that  the  effect  of  the  back  lash  element  in  the  system  dynamics  on 
the  optimal  control  gain  parameters  is  overriding  for  small  load  disturbances. 
In  the  case  that  the  system  is  being  subjected  to  a  load  disturbance  of 
magnitude  somewhat  greater  than  0.05  p.u.  for  the  following  system  parameters 

T  =0.5,  T  =0.5,  M  =  0.04,  G  =  0.01,  and  E  =  0.03 
g  t 

It  was  found  that  the  optimal  control  gain  parameters  with  and  without 
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taking  into  consideration  the  back  lash  effect  are  the  same.  For  load 
disturbances  of  magnitude  somewhat  less  than  0.05  p.u.,  the  optimal 
control  gain  parameters  depend  on  the  disturbance  magnitude.  For  example, 
if  the  system  is  being  subjected  to  a  load  disturbance  of  -  0.004  p.u., 
the  optimal  control  parameters,  when  the  back  lash  element  was  taken  into 
account  in  the  system  dynamics,  were  found  to  be 

k  =  0.0567 
P 

k  =  0.0297 

and  the  cost  functional 

J  =  0.639  x  10~2 

In  the  case  of  applying  the  optimal  control  gain  parameters 

k  =  0.086 
P 

k  =  0.022 

which  are  obtained  by  neglecting  the  nonlinearity  due  to  the  back  lash 
element  (suboptiraal  control  parameters)  to  the  system  with  back  lash,  the 
cost  functional  was  found  to  be 

J5  =  0.725  x  10‘2 


. 


I  - 


149 


Therefore,  by  adopting  a  suboptimal  strategy,  there  is  a  loss  in  the 


optimality  performance  index  as  follows 


e 


(W 


X  100  =  13.45% 


A  comparison  between  the  system  behavior  in  the  two  cases  is  shown 
in  Fig.  (13) . 

To  show  the  effect  of  the  back  lash  element  (governor  dead  band) 
on  the  system  behavior,  the  system  was  considered  to  be  subjected  to 
-  0.01  p.u.,  load  change.  Again,  it  was  found  that  the  optimal  control 
parameters  with  the  back  lash  element  are  different  from  those  obtained 
by  neglecting  it.  Also,  it  was  found  that  the  back  lash  element  introduces 
a  sustained  oscillation  in  system  behavior  of  approximately  0.2  Hz  as  well 
as  increasing  the  maximum  overshoots.  A  comparison  between  the  system 
behavior  in  the  two  cases  is  shown  in  Figs.  (14). 


Conclusion 

The  problem  of  finding  the  optimal  controls  for  a  system  which 
incorporates  a  back  lash  nonlinearity  can  be  handled  by  employing 
Pontryagin's  minimum  principle  and  the  Weierstrass-Erdmann  corner  conditions 
The  problem  of  detecting  the  corner  points  was  solved  by  employing 
interpolation  and  extrapolation  techniques.  The  remaining  part  of  the 
problem  is  to  solve  the  two-point  boundary  value  problem.  This  was  done 
by  employing  the  conjugate-gradient-descent  technique  discussed  in  Chapter  2 


Fig.  (13-a)  Frequency  deviation-time 


Fig.  (13-b)  Backlash  input-time 
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Fig.  (13-c)  Output  power  deviation-time 
Comparison  between  the  system  dynamic  response 
of  steam  area  system  obtained  from  the  nonlinear 
optimal  controller;  optimal  solution, and  the 
linear  optimal  control  parameters;  suboptimal 


solut ion 
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TIME  IN  SECONDS 


Fig.  (14-a)  Frequency  deviation- time 


Fig.  (14-b)  Backlash  input-time 
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Fig.  (14-c)  Output  power  deviation-time 
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TETE  IN  SECONDS 

Fig.  (14-d)  Area  control-t ime 
Comparison  between  the  system  dynamic 
response  of  single  steam  area  with  and 
without  backlash  element;  AL  =  -  0.01  p.u. 


Chapter  6 


Conclusions 

A  general  numerical  method  for  solving  a  deterministic  parameter 
optimization  problem  using  the  steepest  descent-conjugate-gradient 
algorithm  is  presented  in  Chapter  2.  This  general  method  can  be 
implemented  in  large  classes  of  optimal  control  design  problems, 
especially  those  in  which  a  solution  of  a  nonlinear  two-point  boundary 
value  problem  is  required  to  be  found.  A  computer  programme  is  presented 
in  Appendix  A  which  has  been  tested  by  available  solved  dynamic  optimiza¬ 
tion  problems  [65]  . 

The  problem  of  load  frequency  control  (LFC)  of  a  single  and  multiple 
area  interconnected  power  system  was  formulated  as  a  parameter  optimization 
one.  An  especially  attractive  feature  of  the  proposed  scheme  is  that 
it  adopts  present-day  practice  in  the  field  of  load  frequency  control  and 
implements  present —day— opt imal  control  theory.  The  control  law  is  specified 
to  have  proportional  and/or  integral  action  (the  existing  LFC  strategy) . 

By  formulating  the  problem  of  LFC  as  one  of  parameter  optimization  allows 
us  to  compute  the  optimal  control  gain  parameters  off-line.  This  simplifies 
the  on-line  implementation  of  the  proposed  control  strategy. 

In  Chapter  3,  the  problem  of  LFC  of  a  single  area  steam  and  hydro 
power  system  is  considered.  The  optimal  control  parameters  in  the  case 
of  employing  an  integral  and/or  proportional  action  are  found  to  be 
independent  of  the  load  changing.  In  the  case  of  a  single  steam  area 
power  system,  a  proportional-plus-integral-plus-derivative  (PID)  control 
action  is  employed  to  improve  the  system  relative  stability.  Two  schemes 
were  adopted  to  find  the  optimal  control  parameters.  In  the  first  one, 
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the  optimal  control  parameters  were  found  without  any  constraints  on  the 
pole  locations  except  that  they  have  to  be  in  the  L.H.S.  of  the  s-plane. 

In  the  second  one,  the  pole  locations  were  preassigned  in  advance  and 
the  problem  was  converted  to  find  the  optimal  derivative  control  parameter 
(subopt imal  control  parameters).  In  both  cases,  the  optimal  control 
parameters  were  found  to  be  dependent  on  the  load  variations. 

In  spite  of  the  fact  that  the  optimal  or  the  suboptimal  PID  controller 
gives  better  performance  than  the  one  resulting  from  the  optimal  integral 
and/or  proportional  control,  it  is  recommended  to  employ  the  optimal  PI 
controller  to  avoid  the  problem  of  designing  an  adaptive  controller  as 
well  as  the  noise  problem. 

In  Chapter  4,  the  problem  of  LFC  of  an  interconnected  power  system 
is  considered.  First,  the  problem  was  solved  for  the  case  of  neglecting 
the  tie-line  nonlinearity  and  the  areas’  governor  deadband.  In  this  case 
the  optimal  control  parameters  were  found  to  be  independent  of  the  load 
changing.  A  sensitivity  analysis  was  performed,  and  it  was  shown  that 
the  cost  functional  (the  optimality  index)  is  more  sensitive  to  the  areas' 
governor  time  constant  than  the  rest  of  the  areas'  time  constants. 
Therefore,  it  is  recommended  to  put  more  emphasis  on  the  process  of 
identifying  the  areas'  governor-time  constant. 

In  the  case  of  taking  into  account  the  tie-line  nonlinearity,  the 
optimal  control  parameters  were  found  to  be  dependent  on  the  areas' 
load  variations.  However,  the  loss  in  the  system  optimality  when  the 
linear  optimal  control  parameters  are  employed  to  control  the  nonlinear 
system,  is  worthwhile  in  the  sense  that  we  avoid  the  problem  of  designing 
an  adaptive  controller  which  is  required  to  estimate  the  areas'  load 
disturbances  and  accordingly  the  proper  control  parameters. 
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In  Chapter  5,  the  problem  of  taking  into  consideration  a  memory- 
type  nonlinearity  (hysteresis)  in  designing  optimal  control  systems 
is  considered.  The  problem  of  finding  the  optimal  controls  of  a  system 
which  incorporates  a  hysteresis  nonlinearity  can  be  handled  by  employing 
the  Pontryagin’s  minimum  principle  and  the  Weier strass-Erdman  corner 
conditions.  The  problem  of  detecting  the  corner  points  was  solved  by 
employing  extrapolation  and  inverse  interpolation  Hermit e  techniques. 

The  effect  of  the  hysteresis  element  is  prominent  for  small  load  disturbances. 
It  introduces  a  sustained  oscillation  in  the  system  response  and  it  also 
increases  the  maximum  overshoots. 

Suggestions  for  further  study 

The  foregoing  work  can  be  expanded  to  study  the  following: 

(i)  Valve  and  gate  nonlinearities 

In  the  case  of  a  steam  area  system,  the  nonlinearity  due  to  the 
turbine  inlet  valve  arises  from  the  relationship  between  the  governor 
valve  position  x^  and  steam  flow  q  [66] .  As  the  valve  is  progressively 
opened,  the  change  in  steam  flow  for  a  specified  change  in  the  valve 
position  decreases.  There  is  no  available  information  in  the  literature 
about  the  relationship  between  x^  and  q.  Therefore,  one  should  attempt 
to  establish  a  basic  model  for  this  relationship.  Then,  further  study 
could  be  conducted  on  the  problem  of  LFC  taking  into  account  the  valve 
nonlinearity  employing  the  techniques  developed  in  this  thesis. 

In  a  hydro  area  system,  the  relationship  between  the  gate  position 
and  the  turbine  output  power  is  not  linear.  In  the  case  of  a  linear  system 


it  is  given  by 


. 


DP  +  1 
0.5Dp+l 
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where 


P  = 


d_ 

dt 


which  is  a  linear  approximation  of  [67] 


2p  tanh  ptt/2  +  1 
p  tanh  ptt/2  +  1 


in  which 


TTp  =  LV/gH  =  D 


where 


D  =  water  inertia  constant 
V  =  water  speed 
L  =  penstock  length 
H  =  head 

g  =  specific  gravity 

A  similar  relationship  can  be  found  in  [68]. 

(ii)  Boiler  dynamics 

Further  study  can  be  done  by  taking  into  account  boiler  dynamics 
in  the  system  dynamic  equations  [66]. 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 
1  1 
12 

13 

14 

15 

16 

17 

18 

19 

20 
2  1 
22 

23 

24 

25 

26 

27 

28 


46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 


C 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


APPENDIX  A 


This  subroutine  solves  the  problem  of  paramertric  optimization 
using  the  steepest-descent-conjugate-gradient  metods  of 
minimization.  All  the  instructions  for  how  to  use  the  code 
is  documented  in  the  subrotine  grad.  It  is  adivisable  to  test 
the  code  by  one  of  the  problems  which  has  been  discussed  in 
this  desertation. 

SUBROUTINE  GRAD ( YR , FR , YH . FH , TE , YPRED , DFV , SCV , DCV , CV  N1  NC ) 
IMPLICIT  RE AL  *8  (A-H.O-Z) 

REAL *8  YR(N1 ) ,FR(N1 ) ,YH(4,N1 ) , FH( 3 , N 1 ) , TE ( N 1 ) , YPRED(N1 ) , 
A'DFV(NC)  ,  SCV  (  NC  )  .DCV(NC)  ,CV(NC)  ,  F 
INTEGER  RUNGE 

COMMON  Y 1 ( 1 5 1 . 20) . Y2( 1 5 1 , 20) , Y 3 ( 6 . 20 ) . DL 1 , DL 2 
LOGICAL  P 
DATA  P / .  TRUE . / 

EXTERNAL  DE  1  ,DE2 

REFERENCES 

( 1 )  -FOX , RICHARD  "OPTIMIZATION  METHODS  FOR  ENGINEERING  DESIGN" 

ADD  I  SON- WESLEY  COMPANY,  1971 

( 2 )  -D . E . K I RK  "OPTIMAL  CONTROL  THEORY : AN  INTRODUCTION" 

PRENTIC-HALL  INC. ,1970 

(  3) -FLETCHER ,  "PRACTICAL  METHODES  OF  OPTIMIZATION" VOL .  1 
UNCONSTRAINED  OPT IM I ZAT I  ON , JOHN  WILEY  &  SONS,  1980 


29 

C 

SUBROUTINES 

30 

C 

31 

C 

THE  USER  SHOULD  SUPPLY  THE  FOLLOWING  SUBROUTINES 

32 

c 

33 

c 

(  1  )-  DEI  : 

STATE  EQUATIONS 

34 

c 

(2)-  DE2 

CO-STATE  EQUATINOS 

35 

c 

(3)-  SCINT  : 

THE  CO-STATE  EQUATIONS  COFFICIENTS  WHICH  ARE 

36 

c 

FUNCTIONS  OF  THE  STATE  EQUATIONS 

37 

c 

(4)-  T COST  : 

THE  TERMINAL  COST) 

38 

c 

(5)-  COSTAT: 

THE  COSTATE  EQUATIONS  FINAL  CONDITIONS 

39 

c 

40 

c 

N=NO .  OF  THE 

STATE  EQUATIONS 

41 

c 

NC=NO .  OF  THE  CONTROL  PARAMETERS 

42 

c 

N  1  =N+NC+ 1 

43 

c 

NI =MAX IMUM  NUMBER  OF  ITERATIONS 

44 

c 

NT=MAXIMUM  NUMBER  TRIALS  PER  ITERATION 

45 

c 

NR=NUMBER  OF 

ITERATIONS  PERFORMED  BEFORE  RESTSTING  THE 

CONJUGATE-GRADIENT  METHOD ( 1  GIVES  THE  STEEPEST-DESCENT ) 
RAT  1 0  =  M INI MUM  VALUE  OF  ABS ( D  H( ALPHA ) /D  H(BETA)) 
DFMIN=MINIMUM  VALUE  OF  ABS ( THE  GRADIENT  OF  COST  FUNCTIONAL 
W.R.T.  THE  CONTROL  VECTOR) 

FMIN=LOWER  BOUND  ON  THE  COST  FUNCTIONAL 
FMAX=UPPER  BOUND  ON  THE  COST  FUNCTIONAL 
ALMIN=MINIMUM  INITIAL  VALUE  OF  ALPHA 
ALMAX  =  MAX I  MUM  INITIAL  VALUE  OF  ALPHA 
RMIN=MINI MUM  FRACTION  OF  (ALPHA-BETA) 

REDUCT  I  ON ( < 1 ) 

RMAX=MAXIMUM  FRACTION  OF  (ALPHA-BETA)  USED  IN 
EXTERAPOLATION(> 1 ) 

XMAX=F INAL  TIME 
H= INTEGRATION  STEP  SIZE 
F  = J ( C )  : THE  VALUE  OF  THE  COST  FUNCTION 


USED  IN  INTERVAL 
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6  1 

62 

63 

64 

65 

66 

67 

68 

69 

70 

7  1 

72 

73 

74 

75 

76 

77 

78 

79 

80 

8  1 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

1  10 

1  1  1 

1  12 

1  13 

1  14 

1  15 

1  16 

1  17 

1  18 

1  19 

120 


160 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DF  V 1 NC ) =  DU ( C  )  /DC 

ARRANGE  THE  EQUATIONS  IN  DEI  AS  FOLLOWS: 
-THE  STATE  DYNAMIC  EQUATIONS 
-THE  CONSTRAIN  DYNAMIC  EQUATIONS 
-THE  CONTROL  DYNAMIC  EQUATIONS 
-THE  COST  FUNCTIONAL  DYNAMIC  EQUATION 

ARRANGE  THE  EQUATIONS  IN  DE2  AS  FOLLOWS: 
-THE  CO-STATE  DYNAMIC  EQUATIONS 
-THE  CO-CONSTRAIN  DYNAMIC  EQUATIONS 
-THE  CO-CONTROL  DYNAMIC  EQUATIONS 
-THE  CO-COST  FUNCTIONAL  DYNAMIC  EQUATION 


10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
1  10 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


( 


i  H '  ) 


1 

120  FORMAT 


1 


130 

140 

150 

160 


FORMAT 

FORMAT 

FORMAT 

FORMAT 


XMAX 

( 2D  1 8 .  10) 

(4D  18 .  10) 

(  4X  ,  4D18 . 10) 

(41 10) 

(  '  1  '  ) 

'  NI , NT , NR , NC ' ) 

RATIO .DFMIN , FMIN, FMAX  '  ) 

ALMIN, ALMAX, RMIN, RMAX' ) 

THE  CONTROL  VECTOR ' ) 

'  ITERATION^',  12,  ' 

10) 

ALPHA  = '  ,  D18 .  10.  '  DF  = 

10) 

'  RESET  THE  CONTROL  VECTOR 

'  DESIGN  COMPLETE ' ) 

'  ITERATION  LIMIT ' ) 

'  TRIAL  LIMIT ' ) 


(//. 

(' 

(  ' 

(  ' 

(//. 
D  1  8  . 
(  ' 

D  18  . 
(//. 
(//. 
(//. 
(//, 


TRI AL= 


12 


D 1 8 . 10, 


WRITE  (6,60) 

READ  (5,20)  XMAX.  H 
(5,50)  NI,  NT. 


READ 

READ  (5,30)  RATIO, 
READ  (5,30)  ALMIN. 
WRITE  (6,10) 

WRITE  (6.20)  XMAX,  H 
WRITE  (6,70) 

WRITE  (6,50) 

WRITE  (6,80) 

WRITE  (6,40) 

WRITE  (6,90) 

WRITE  (6.40) 


NR,  NC 
DFMIN, 
ALMAX, 


FMIN, 
RMIN , 


FMAX 

RMAX 


NI,  NT,  NR,  NC 


RATIO,  DFMIN,  FMIN,  FMAX 


ALMIN.  ALMAX,  RMIN.  RMAX 


N  =  NI  -  NC 
NIC  =  O 
NT  C  =  O 
NRC  =  NR 


-  1 


C  - 

C  DC=( ALPHA-BETA ) 

C  - 

C  INITIALIZE  THE  CONTROL  VECTOR 

C  - 

DO  170  J  =  1 .  NC 
170  CV( J )  =  Y 1 ( 1 . J  +  N) 

C 

C  =  0 . DO 


F  =  '  , 
1  DF V 1  =  '  , 


12  1 

122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

14  1 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

17  1 

172 

173 

174 

175 

176 

177 

178 

179 

180 
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X  =  1 . DO- 1 . DO  /  RMAX 

C  - 

c  INTEGRATE  THE  STATE  EQUATIONS  AND  THE  COST  FUNCTION 

C  DYNAMINC  EQUATIONS  FORWARD 

C  - 

XI  =  XMAX 
X2  =  0 . DO 

II  =  (XMAX  +  .  5  1  H  )  /  H  +  1 

CALL  HAM I  MG ( DEI,  YR,  FR ,  YH.  Y1,  FH,  TE,  YPRED ,  N1,  XI, 


C  COMPUTE  THE  TERMINAL  COST 

C  - 

DO  ISO  J  =  1 ,  N 1 
180  Y  F!  (  d )  =  Y  1(  I  1  .  J) 

CALL  T  COST ( YR ,  FT) 

F  =  YI(II.NI)  +  FT 


C  INTEGRATE  THE  CO-STATE  EQUATIONS  BACKWARD 

C  COMPUTE  THE  CO-STATE  TERMINAL  CONDITIONS 


CALL  COSTAT ( YR ) 

DO  190  d  =  1 ,  N1 
190  Y2( 1 ,d)  =  YR ( J ) 


XI  =  XMAX 

C  - 

C  GENERATE  THE  SOLUTIONS  FOR  THE  STATE  EQUATIONS  AT ( XMAX , 

C  XMAX-H/2 . XMAX-H . XMAX-5H/2 ) ; US ING  RUNGE-KUTTA. 

c  - 


c 

c 

c 

c 

c 


c 


DO  200  J  =  1 ,  N 1 
Y3(  1  , J)  =  Y 1 ( I  1  . J) 

200  YR ( d )  =  Y 1 ( I  1 , d ) 

I  =  1 

210  K  =  RUNGE(N1 , YR . FR, XI , -H/2 . ) 
IF  (K  .NE.  1)  GO  TO  220 
CALL  DE1(YR,  FR.  XI) 

GO  TO  210 
220  1=1+1 

DO  230  d  =  1 ,  N1 
230  Y3 ( I , d )  =  YR ( d ) 

IF  (I  .LT.  6)  GO  TO  210 
XI  =  0 . DO 
X2  =  XMAX 
P  =  .TRUE. 


CALL  HAMING( DE2 . 


YR.  FR. 


YH,  Y  2 . 


FH,  TE.  YPRED,  N1, 


COMPUTE  THE  GRADIENT  OF  THE  COST  FUNCTIONAL  W.R.T. 
CONTROL  VECTOR. 


DO  240  d  =  1 ,  NC 
240  DFV(d)  =  Y2 ( I  1 , N  +  d) 


XI,  - 


DF  =  O . DO 
DO  250  d  =  1 .  NC 
SCV(d)  =  CV(d) 

DCV(d)  =  -DFV(d) 

250  DF  =  DF  +  DCV(d)  *  DFV(d) 
C 


,  X2 ,  P) 


H,  X2 ,  P) 


-as 
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18  1 
182 

183 

184 

185 

186 

187 

188 

189 

190 
19  1 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 
209 
2  10 
2  1  1 
2  12 
2  13 
2  14 
215 


C 

C 


DF  =  -<DF ( C  )/DC . DF( C )/DC> 
SVP  =  -DF 

SDF  =  RATIO  *  DABS ( DF ) 
D3  =  DSQRT(SVP) 

WRITE  (6,110)  NIC,  NTC 
WRITE  (6. 120) 

WRITE  (6.100) 

WRITE  (6,40)  ( CV( I ) 

IF  ( FMAX  . GT.  F)  GO 
WRITE  (6, 130) 

RETURN 


C.  DF,  D3 


I = 1 . NC ) 
TO  260 


260 


CHECK  THE  TERMINATION 

03 )  GO  TO  270 


OF  THE  TRIALS  AND  THE  ITERATINS 


IF  (DFMIN  ,LT 
WRITE  (6, 140) 

WRITE  (6,60) 

WRITE  (6,30)  ( CV( I ) 
RETURN 


1= 1 ,NC) 


C 

c 

c 

c 


INTIALIZE  THE 
ONE -D I  MENS  I OL 


INCREMENT  OF 
DIFFERENTIAL 


THE  STEP  SIZE  DC 
SEARCH  ALGORITHM. 


USING  THE 


270 

280 

290 


2. DO  *  ( FMIN  -  F)  /  DF 
DM IN  1 ( ALMIN , DM AX  1 ( DC , ALMAX ) ) 


300 


C 

C 


DC  = 

DC  = 

NTC  =  NTC  +  1 
A  =  C 
FA  =  F 
DF  A  =  DF 
C  =  A  +  DC 
DO  300  I  =  1  ,  NC 
CV( I )  =  SCV(I)  + 


DCV ( I ) 


INTEGRATE  THE 


2  16 

C 

CONTROL  VECTOR 

217 

C 

-  - - 

2  18 

DO  310  J  =  1 ,  NC 

219 

310 

Y  1  (  1  , d  +  N)  =  CV( J) 

220 

XI  =  XMAX 

22  1 

X  2  =  0 . DO 

222 

P  =  TRUE. 

223 

CALL  HAMING(DE 1 ,  YR 

224 

DO  320  J  =  1 ,  N 1 

225 

320 

YR ( J )  =  Y  1  ( I  1  , J) 

226 

CALL  TCOST ( YR ,  FT) 

227 

F  =  Y1(I1,N1)  +  FT 

228 

CALL  COST  AT ( YR ) 

229 

DO  330  <J  =  1  ,  N 1 

230 

330 

Y2(  1,J)  =  YR ( U ) 

231 

XI  =  XMAX 

232 

DO  340  J  =  1 ,  N 1 

233 

Y3(  1  . J)  =  Y1  ( I 1 ,J) 

234 

340 

YR ( d )  =  Y 1 ( I  1  . J) 

235 

I  =  1 

236 

350 

K  =  RUNGE(N1 . YR.FR.X 

237 

IF  (K  . NE.  1  )  GO  TO 

238 

CALL  DE 1 ( YR ,  FR,  XI) 

239 

GO  TO  350 

240 

360 

1  =  1  +  1 

AND  THE  CO-STATE  EQUATIONS  WITH  THE 


NEW 


FR,  YH,  Y 1 ,  FH.  TE.  YPRED ,  N1,  XI.  H,  X2,  P) 


1 , -H/2 
360 


. 


. 


241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 

261 

262 

263 

264 

265 

266 

267 

268 

269 

270 

27  1 

272 

273 

274 

275 

276 

277 

278 

279 

280 

281 

282 

283 

284 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 


n  o  o  o  n  o 
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c 

c 

c 


DO  370  J  =  1 ,  N1 
370  Y3( I . J)  =  YR( J) 

IF  ( I  .LT.  6)  GO  TO  350 
XI  =  O . DO 
X2  =  XMAX 
P  =  .TRUE. 

CALL  HAMI NG( DE  2 .  YR,  FR,  YH,  Y2.  FH,  TE .  YPRED .  N1  XI 
DO  380  U  =  1  ,  NC 
380  DFV(U)  =  Y2( I  1 , N  +  J) 


-H.  X2 , 


CONTROL  VECTOR (CV=CV+ALPHA*GRAD.  OF  THE  COST  FUNCTIO  W.R.T  CV) 


DF  =  0 . DO 
D 1  =  0 . DO 
DO  390  J  =  1 ,  NC 

DF  =  DF  +  DCV(U)  *  DFV(U) 
390  D 1  =  D 1  +  DFV(U)  *  DFV(J) 

D2  =  DSQRT(DI) 


CHECK  FOR  THE  NORMAL  TERMINATION  OF  THE  TRIAL  PROCEDURE: 
<DIRECTION( I ) , -Dd(CV( I ) ) >/<DIRECTION( I ) , -Dd(CV( 1+1 ) )> .LE .RATIO 


IF  ( SDF  .GE.  DABS ( DF ) )  GO  TO  630 
WRITE  (6.1 10)  NIC,  NT  C ,  F 
WRITE  (6. 120)  C,  DF.  D2 
WRITE  (6.100) 

WRITE  (6,40)  (CV( I) , 1=1 ,NC) 

IF  ( FMAX  . GT.  F)  GO  TO  420 
400  IF  (NT  .GT.  NT C )  GO  TO  4  10 
GO  TO  630 


C  - 

C  WHENEVER  F ( CV+STEP*D I RCET ION( I ) ) . GT . FMAX  REDUCE  STEP  BY  A 

C  CERTAIN  FACTOR  RMIN  IN  THE  RANGE  OF  (0,1)  AND  ADDITOINAL  TRIAL 

C  WILL  BE  PERFORMED 

C  - 


410  DC  =  RMIN  *  DC 
C  =  A 
F  =  FA 
DF  =  DF  A 
GO  TO  290 


P) 


CHECK  FOR  <DIRECTI ON ( I ) , F (CV+STEP*DIRECTION( I )> . GT .0 


420  IF  (DFMIN  .LT.  D2 )  GO  TO  430 
WRITE  (6.140) 

GO  TO  660 

430  IF  (DF)  440.  440.  470 
440  IF  (NT  . LT.  NT C )  GO  TO  630 


ESTIMATE  THE  OPTIMUM  STEP  SIZE  BY  QUADRATIC  EXTERAPOLATION 


IF  ( X +  D F A  .LE.  DF)  GO  TO  450 
C  DC=DC+RMAX 

C  GO  TO  1021 

C  - 

C  BY  BASS  THE  QUADRATIC  EXTRAPOLATION  AND  START  NEW  TRIAL 

C  - 

GO  TO  630 

450  DC  =  DC  ‘  DF A  /  ( DFA  -  DF) 

GO  TO  290 


' 


-• 


301 

302 

303 

304 

305 

306 

307 

308 

309 
3  10 
3  1  1 
3  12 
313 
3  14 
315 
3  16 
317 
3  18 

319 

320 
32  1 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 
34  1 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 


460  WRITE  (6, 160) 
GO  TO  660 


C 

c 


480 


C 

c 

c 


ESTIMATE  TIHE  OPTIMUM  STEP  SIZE  BY  QUBIC  INTTERAPOLATION 

470  B  =  C  ~~ 

FB  =  F 
DFB  =  DF 

01  =  DF A  +  DFB  +  3.00  *  ( FA  -  FB)  /  DC 
D2  =  2 . DO  *  ( DF  A  +  D 1 ) 

D3  =  DFA  +  DFB  +  2  DO  *  D 1 
IF  (1.D-5*DABS(D2)  .LE.  DABS ( D3 ) )  GO  TO  490 
C  -  B  -  (DFA  +  2 . DO*D 1 )  *  DC  /  D2 
GO  TO  500 

D2  =  DSQRT ( D  1  **2  -  DFA+DFB) 

C  =  B  -  (DFB  +  D1  -  D2 )  *  DC  /  D3 
DO  510  I  =  1 ,  NC 
CV(I)  =  SCV(I)  +  c  *  DCV(I) 

INTEGRATE  THE  STATE  AND  CO-STATE  EQUATIONS 


490 

500 
5  10 


520 


NC 

CV(J) 


FR,  YH,  Y 1 ,  FH.  TE .  YPRED.  N1.  XI,  F 


530 


540 

550 

560 

570 

580 


590 


DO  520  d  =  1 
Y  1  (  1 , d  +  N)  ■ 

XI  =  XMAX 
X2  =  0 . DO 
P  =  .TRUE. 

CALL  HAMING( DEI,  YR 
DO  530  d  =  1 ,  N 1 
YR  (  d )  -  Y 1 ( I  1 . d ) 

CALL  T  COST ( YR ,  FT) 

F  =  Y  1  ( I  1 , N 1 )  +  FT 
CALL  COSTAT ( YR ) 

DO  540  d  =  1 ,  N 1 
Y2( 1 ,d)  =  YR ( d ) 

XI  =  XMAX 
DO  550  d  =  1 ,  N 1 
Y 3 (  1  ,  d )  =  Y 1 ( I  1 , d) 

YR ( d )  =  Y 1 ( I  1 , d) 

I  =  1 

K  =  RUNGE(N1 , YR , FR, XI . -H/2 . ) 

IF  (K  .NE .  1 )  GO  TO  570 

CALL  DE  1  ( YR ,  FR.  XI) 

GO  TO  560 
1  =  1  +  1 
DO  580  d  =  1 ,  N1 
Y3( I , d )  =  YR ( d ) 

IF  (I  .LT.  6)  GO  TO  560 
XI  =  O . DO 
X2  =  XMAX 
P  =  .TRUE. 

CALL  HAMING( DE2 ,  YR.  FR.  YH,  Y2.  FH,  TE.  YPRED  N1  XI 
DO  590  d  =  1,  NC  * 

DFV(d)  =  Y2( I  1 ,N  +  d) 


DF  =  O . DO 
D  1  =  0 . DO 
DO  600  I  =  1 ,  NC 

DF  =  DF  +  DCV ( I )  *  DF V ( I ) 
600  D 1  =  D 1  +  DFV( I )  +  DFV(I) 

D2  =  DSQRT ( D 1 ) 


X2 ,  P) 


X2 ,  P) 


>K  .  t  -  i  <  U  » 


36  1 

362 

363 

364 

365 

366 

367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 

403 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1  1 

12 

13 

14 

15 

16 

17 


non 
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IF  ( SDF  .GE.  DABS ( DF ) )  GO  TO  630 
NT  C  =  NT  C  +  1 
WRITE  (6,110)  NIC ,  NT  C ,  F 
WRITE  (6, 120)  C,  DF,  D2 
WRITE  (6,100) 

WRITE  (6,40)  (CV(  I  )  . 1  =  1  ,NC) 

IF  ( FMAX  .LE.  F)  GO  TO  400 
IF  (NT  .LE.  NT C )  GO  TO  630 
IF  ( DF )  620,  610.  610 
610  DC  =  C  -  A 
GO  TO  470 
620  A  =  C 
FA  =  F 
DF A  =  DF 
DC  =  B  -  C 
GO  TO  480 
630  NIC  =  NIC  +  1 
DC  =  C 
C  =  O.ODO 
NT  C  =  O 

IF  ( NRC  .GT.  NIC)  GO  TO  640 


RESTART  THE  CONJUGATE  GRADINET  PROCEDURE 


SVP  =  1 . D50 
NRC  =  NRC  +  NR 
640  B  =  D 1  /  SVP 
SVP  =  D 1 
DF  =  0 . DO 
DO  650  I  =  1 .  NC 
SCV(I)  =  CV(I) 

DCV(I)  =  B  *  DCV(I)  -  DFV(I) 

650  DF  =  DF  +  DCV(I)  *  DF V( I ) 

SDF  =  RATIO  *  DABS ( DF ) 

WRITE  (6,  1  10)  NIC,  NT  C ,  F 
WRITE  (6, 120)  C,  DF .  D2 
WRITE  (6. 100) 

WRITE  (6,40)  (CV( I ) , 1=1 , NC ) 

IF  ( N I  .GT.  NIC)  GO  TO  280 
WRITE  (6. 150) 

660  WRITE  (6,60)  ( SCV ( I ) , I = 1 , NC ) 

RETURN 
END 

FUNCTION  RUNGE ( N ,  Y,  F.  X.  H) 

C  The  FUNCTION  RUNGE  EMPLOYS  THE  FOURTH-ORDER  RUNGE-KUTTA 

C  METHOD  WITH  CUTTA ' S  COEFFICIENTS  TO  INTEGRAT  A  SYSTEM  OF 

C  N-SIMULTANEOUS  FIRST  ORDER  DI FFERENTIONAL  EQUATIONS. 

C  Maximun  number  of  equations  is  equal  to  50 

C 

IMPLICIT  REAL*8(A  -  H,0  -  Z) 

REAL *8  X.  H 
INTEGER  RUNGE 

RE AL *8  PHI (50) ,  SAVEY(50),  Y(50).  F(50) 

DATA  M  /O/ 

C 

M  =  M  +  1 

GO  TO  ( lO.  20.  40.  60.  80).  M 
C 

C  PASS  1 


18 

19 

20 

2  1 

22 
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32 
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37 
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42 

43 
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50 

51 
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55 

56 
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58 
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62 
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c  - 

c 

10  RUNGE  =  1 
RETURN 
C 

c  PASS  2 
C  - 

c 

20  DO  30  J  =  1 ,  N 

SAVE  Y ( U )  =  Y(U) 

PHI(J)  =  F ( J ) 

30  Y(J)  =  SAVE Y ( J )  +  5  ♦  H  *  F(J) 

X  =  X  +  .5  *  H 
RUNGE  =  1 
RETURN 
C 

C  PASS  3 

C  - 

c 

40  DO  50  J  =  1 ,  N 

PHI(U)  =  PHI ( J )  +  2.  *  F  (  J ) 

50  Y ( U )  =  SAVE Y ( U )  +  . 5  *  H  *  F(U) 
RUNGE  =  1 
RETURN 

PASS  4 


60  DO  70  U  =  1 .  N 

PHI ( J )  =  PHI(J)  +  2.  *  F  (  J ) 

70  Y(U)  =  SAVE Y ( U )  +  H  *  F(J) 

X  =  X  +  .5  +  H 
RUNGE  =  1 
RETURN 
C 

C  PASS  5 

C  - 

c 

80  DO  90  J  =  1  ,  N 

90  Y(J)  =  SAVE Y ( J )  +  (PHI(J)  +  F(J))  *  H  /  6 . 

M  =  O 
RUNGE  =  0 
RETURN 
C 

END 

C 

FUNCTION  HMING( N ,  Y,  YPRED ,  F,  X,  H,  TE,  PRED) 
c  - 

C  HMING  IMPLEMENTS  HAMMING'S  PREDICTOR-CORRECTOR  ALGORITHM 

C  TO  SOLVE  N  S IMLUT ANEOUS  FIRST-ORDER  DI F FR ENT I ONAL EQUATIONS . 

C 

IMPLICIT  RE AL K  8 ( A  -  H.O  -  Z) 

REAL  +8  X,  H 
INTEGER  HMING 
LOGICAL  PRED 

REAL  *8  YPRED ( N ) ,  TE ( N ) ,  Y(4.N).  F(3,N) 

C 

C  is  CALL  FOR  PREDICTOR  OR  CORRECTOR  SECTION 

c  - 

c 


■ 
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82 
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99 
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114 
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120 
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IF  (  .NOT.  PRED)  GO  TO  40 
PREDICTOR  SECTION  OF  HMING 

COMPUTE  PREDICTED  Y(d)  VALUES  AT  NEXT  STEP 


C 

C 

c 


c 

c 

C 

c 

c 


c 

c 

c 

c 


DO  10  J  =  1  ,  N 

10  YPRED(U)  =  Y ( 4 , J )  +  4  .  *  H  *  ( 2  .  *  F (  1  U) 
UPDATE  THE  Y  AND  F  TABLES 


F(2.J)  +  2.*F(3,J))  /  3. 


DO  20  J  =  1 .  N 
DO  20  K  =  1 ,  3 
L  =  5  -  K 

Y(L.U)  =  Y(L  -  1.J) 

20  IF  (L  .LT.  4)  F(L.d)  =  F ( L  -  1  j) 


Y(J)  VALUES  USING  THE  TRUNCATION  ERROR 
ESTIMATES  FROM  PREVIOUS  STEP,  INCREMENT  X  VALUE 


DO  30  d  =  1 .  N 

30  Y  (  1  ,  d )  =  YPRED(U)  *-  112.  *  TE(d)  /  9 
X  =  X  +  H 


SET  PREDAND  REQUEST  UPDATE  DERIVATIVE  VALUES 


C 

c 

c 

c 

c 

c 


PRED  =  .FALSE. 
HMING  =  1 
RETURN 


CORRECTOR  SECTION  OF  HMING 
COMPUTE  CORRECTED  AND  IMPROVED 
TRUNCATION  ERROR  ESTIMATES  FOR 


VALUES  OF  THE  Y(d)  AND  SAVE 
THE  CURRENT  STEP 


40  DO  50  d  =  1 .  N 

YO.d)  =  (9.*Y(2,d)  -  Y  (  4  ,  d )  +  3  .  *H*  (  F  (  1  ,  d )  +  2.*F(2,d)  -  F(3,d) 
1  ))/  8 . 

TE(d)  =  9.  +  ( Y ( 1 , d )  -  YPRED(J))  /  121. 

50  Y(1,d)  =  Y (  1  , d )  -  TE  (  J ) 


SET  PRED  AND  RETURN  WITH  THE  SOLUTION  FOR  THE  CURRENT  STEP 


C 

C 


c 


PRED  =  .TRUE. 
HMING  =  2 
RETURN 

END 


SUBROUTINE  HAMING( DE ,  YR, 
1  P) 


FR.  YH,  Y,  F,  TE,  YPRED ,  N,  XMAX ,  H,  X, 


IMPLICIT  REAL*8(A  -  H.O  -  Z) 

INTEGER  RUNGE.  HMING 
REAL* 8  X,  H.  XMAX 

RF  AL  *8  TE ( N) ,  YR ( N ) ,  FR(N),  Y(101,N),  F(3,N),  YPRED(N),  YH(4  N) 
LOGICAL  P 
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139 
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150 
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160 
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163 
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165 
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168 
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171 
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173 
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177 

178 
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c 

c 

c 

c 

c 

c 

c 


external  de 


Ma x 1 "  no  of  output  points" is  101(11 
any  number  by  changing  the  dimension 


could  changed  into 

of  matrix  Y( * , N) ) 


SET  INITIAL  TRUNCATION  ERROR 
SET  THE  FIRST  ROW  OF  V  MATRIX 


TO  ZERO 


I  =  1 

DO  10  d  =  1  N 
TE(J)  =  O. 

YH(4,d)  =  Y(I,J) 

10  YR(J)  =  Y(I,J) 


CALL  RUNGE  TO  INTEGERATE  ACROSS  THE  FIRST  THREE  STEPS 

20  K 1  =  RUNGE (N.YR.FR.X.H) 

IF  ( K 1  . NE .  1 )  GO  TO ’ 30 

CALL  DE ( YR ,  FR.  X) 

GO  TO  20 
30  I  =  I  +  1 

DO  40  U  =  1 ,  N 
K  =  5  -  I 
YH(K.J)  =  YR(  J) 

40  Y( I , J)  =  YR(J) 

CALL  DE ( YR .  FR,  X) 

DO  50  I  1  =  1 ,  N 
50  F ( K ,  I  1  )  =  FR(  I  1  ) 

IF  (I  . LT.  4)  GO  TO  20 

CALL  HMING  TO  INTEGEATE  ACROSS  THE  NEXT  STEPS 

60  M  =  HMING(N. YH. YPRED.F  X  H  TE  P) 

DO  70  12  =  1.  N 
70  YR( 12)  =  YH( 1,12) 

CALL  DE ( YR ,  FR,  X) 

DO  80  13  =  1,  N 
80  F(  1  . 13)  =  FR( 13 ) 

XMAX)  . L E .  H/2.)  GO  TO  100 
XMAX)  . GE .  H/2.)  GO  TO  100 


DO  90  U  =  1 ,  N 
90  Y(I.J)  =  YR(J) 
GO  TO  60 

100  RETURN 
END 


IF  (H  .LT.  O . DO  .AND.  (X  - 
IE  (H  .GT.  0 . DO  .AND.  (X  - 
IF  (M  . EO.  1  )  GO  TO  60 
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32 

33 
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non 


APPENDIX  B 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


THIS  SUBROUTINE  SIMULATE  THE  BACKLASH  ELEMENT  DIGITIALLY 

SUBROUTINE  BKL ASH( YBL I ,  YPBLI.  YBLO ,  X,  H,  DB  YI  YIP  YS  YSP 
1  N.  M,  STATUS)  .  x  ,  is,  t5k. 

IMPLICIT  REAL*8( A  -  H.O  -  Z) 

INTEGERYSTATUSP(2)  ’  YR(2)’  YI(N)’  YIp(N).  YS (N ) ,  YSP(N) 

COMMON  /BL0CK3/  DUMMY,  M0DE(2),  MODEDB 
COMMON  /BL0CK4/  K(301) 

COMMON  /BL0CK5/  U1 1 
LOGICAL  START 


Y ( * )  WORK  VECTOR  OF  LENGTH  2  TO  STORE  2-VALUES  OF  THE  INPUT 
TO  THE  BACKLASH  BEFORE  REVERSING  THE  DIRECTION 
YP ( * )  WORK  VECTOR  OF  LENGTH  2  TO  THE  DERIVATIVE  OF  Y(*) 

WITH  RESPECT  TO  TIME 

YI(»)  WORK  VECTOR  OF  LENGTH  N  TO  STORE  THE  VALUES  OF  THE 
PREVAILING  SYSTEM  STATES  AT  X 
Y I P ( * )  WORK  VECTOR  OF  LENGTH  N  TO  STORE  THE  VALUES  OF  THE 
PREVALING  SYSTEM  DERIVALTIVES  AT  X 
YS (  * )  WORK  VECTOR  OF  LENGTH  N  TO  STORE  THE  VALUES  OF  THE 
SYSTEM  STATES  AT  X-H 

YSP ( * )  WORK  VECTOR  TO  STORE  THE  DERIVATIVE  OF  YS(*) 

YBLI  INPUT  TO  THE  BACK-LASH 

YPBLI  THE  DERIVATIVE  OF  THE  INPUT  TO  THE  BACK-LASH 
YBLO  OUTPUT  FROM  THE  BACK-LASH 

N  NUMBER  OF  STATE  EQUATIONS 

M  THE  INDEX  OF  INPUT  STATE  TO  THE  BACKLASH  ELEMENT 

K ( * )  VECTOR  OF  DIMENSION  IMAX 

K l * ) =0  DEAD  BAND  IS  ACTIVE 
K (  *  )  =  1  DEAD  BAND  IS  INACTIVE 
MODE ( 2 ) :  PRESENT  MODE  OF  OPERATION 
MODEM):  PAST  MODE  OF  OPERT A 1  ON 
+1  RISING  MODE 
-  1  FALLING  MODE 
MODEDB  : 1  DEAD  BAND  IS  ACTIVE 

: 0  DEAD  BAND  IS  INACTIVE 
STATUS  : 2  INTERPOLATION  OR  EXTRAPOLATION  MODE 
: 1  OTHERWISE 
START  = . TRUE  . 

MODE  OF  OPERATION  IS  AT  THE  HALF  OF  DEAD-BAND 
(INITIAL  CONDITION) 

START  = . FALSE  . 

MODE  OF  OPERATION  WITH  FULL  DEAD-BAND 
YR ( * )  THE  LAST  TWO  POINTS  OF  REVERSING  DIRECTION 
DB  DEAD  BAND 

J  1  -ST ATUS  OF  THE  PREVAILING  TWO  POINTS 


INITIALIZE  THE  SYSTEM  MODE  OF  OPERATION 


STATUS  =  1 

IF  (X  . GT .  1 . 1*H)  GO  TO  20 

START  =  .TRUE. 

MODEDB  =  1 
XR  =  0.0 


169 


61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

7  1 

72 

73 

74 

75 

76 

77 

78 

79 

80 

8  1 

82 
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92 
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99 
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1  1  1 

112 

1  13 

1  14 

115 

1  16 

1  17 

1  18 

1  19 

120 
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C  INITIALIZE  Y.YP.YR.K 

C  - 

DO  10  I  =  1 .  2 
MODE (  I  )  =  0 
Y ( I )  -  0.0 
YP ( I )  =  0.0 
10  YR ( I  )  =  0.0 
K  (  1  )  =  O 

c 


c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


GO  TO  40 


CHECK  UP  FOR  THE  MODE  OF  OPERATION 


20  IF  (YPELI)  70.  30,  110 
30  IF  (YELI  .LT.  0  .AND.  M0DE(2) 
IF  ( YBL I  .LT.  0  .AND.  M0DE(2) 
IF  (YELI  .GT.  0  .AND.  MODE (2) 
IF  (YBL I  .GT.  0  .AND.  M0DE(2) 
MODEDE  =  1 
GO  TO  220 


LT  . 

0) 

YBLO  = 

YBL I  -  DB 

/ 

2 

GT  . 

0) 

YBLO  = 

YBL I  +  DB 

/ 

2 

LT  . 

0) 

YBLO  *= 

YBL I  -  DB 

/ 

2 

GT  . 

0) 

YBLO  = 

YBL I  +  DB 

/ 

2 

AT  THE  FIRST  STEP  OF  INTEGERAT ION  (X.EQ.H)  THE  MODE 
ASSME  THAT  AT  X=H  THAT  YBL0=0.0  OR  YBLI+-BD/2. 


40  IF  (YPELI  .  LT  .  0)  M0DE(2)  =  -1 
IF  (YPBLI  .GT.  0)  MODE ( 2 )  =  1 
IF  (DABS(YBLI)  .LE.  DB/2.)  GO  TO  50 
IF  ( YBL I  .LT.  0)  YBLO  =  YBL I  +  DB  /  2 . 
IF  (YELI  .GT.  0)  YBLO  =  YBL I  -  DE  /  2 . 
MODEDB  =  0 
START  =  .FALSE. 

GO  TO  60 
50  YBLO  =  0.0 
60  GO  TO  220 


FALLING  MODE  OF  OPERATION 


70  MODE ( 1 )  =  MODE (2 ) 

MODE ( 2  )  =  -1 

<J1  =  MODE  (  1  )  *  MODE  (  2  ) 

IF  (MODEDB  .EO.  1)  GO  TO  180 
IF  ( J1  .LT .  0)  GO  TO  140 
IF  (  .NOT.  START)  GO  TO  90 
IF  (DAES(YELI)  -  DB/2.)  80.  80,  90 
80  YBLO  =  O.ODO 
GO  TO  220 

90  IF  (YELI  .GT.  O)  YBLO  =  YBL I  -DB/2. 
IF  ( YBL  I  LT.  0)  YBLO  =  YBL I  +  DB  /  2  . 

100  FORMAT  (2X,  E14.5.  213,  4E14.5/) 

GO  TO  220 


C  - - - 

C  RISING  MODE  OF  OPERATION 

L  - 

1 10  MODE ( 1 )  =  MODE ( 2 ) 

MODE (2)  =  1 

J1  =  MODE ( 1 )  *  MODE (2) 

IF  (MODEDB  .EO.  1)  GO  TO  180 

IF  ( U 1  . LT .  1 )  GO  TO  140 

IF  (  .NOT.  START)  GO  TO  130 

IF  ( DABS ( YBL I  )  -  DB/2.)  120,  120,  130 


' 
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142 

143 
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15  1 

152 
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154 
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156 

157 
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159 

160 
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120  YBL0  =  0.0D0 
GO  TO  220 

130  IF  (YBLI  .LT.  0)  YBLO  =  YBL I  ♦DB/2. 
IF  (YBLI  .GT.  0)  YELO  =  YBLI  -DB/2 
GO  TO  220 
C 

140  YR  (  1  )  =  YR ( 2 ) 


IMPEMENT  THE  THIRD-ORDER  HERMITE  INTERPOLATION  FORMULA 
TO  FIND  THE  REVERSING  POINTS 


CALL  I NTR ( Y ,  YP .  YS,  YSP.  X,  H,  YR(2),  XDB ,  XR,  N.  M) 
STATUS  =2 


CHECK  FOR  THE  MODE  OF  OPERATION 


C 


c 


c 


IF  ( DABS ( YS ( M )  -  YR(2))  .LT.  DB )  MODEDB  =  1 
IF  ( DABS ( YS ( M  )  -  YR ( 2  )  )  .GE.  DB  )  MODEDB  =  0 

YP( 1 )  =  YP ( 2  ) 

Y  (  1  )  =  Y  (  2  ) 

YP ( 2 )  =  YSP(M) 

Y  (  2  )  =  Y  S  (  2  ) 


IF  (  .NOT.  START)  GO  TO  150 
IF  (START  .AND. 

150  IF  ( MODE ( 1 )  .GT . 

IF  ( MODE ( 1 )  .GT. 

IF  ( MODE ( 1 )  .LT. 

IF  ( MODE (  1  )  .LT. 

MODE ( 1 )  =  MODE ( 2 
IF  ( YP( 2  )  .  LT  .  0)  MODE ( 2  )  =  -  1 
IF  ( VP ( 2 )  GT.  0)  MODE ( 2  )  =  1 
IF  ( DABS ( YR ( 2  )  -  V  R (  1 ) )  .LT.  DB )  GO  TO  250 


DABS ( YR ( 2 )  -  YR (  1  ) ) 

.LE 

.  DB/2. 

)  GO  TO  160 

0 

.AND  . 

YR  (  2 ) 

.GT  . 

0.  ) 

YBLO  = 

YR ( 2 )  -  DB 

/ 

2 

0 

.  AND  . 

YR  ( 2  ) 

.LT 

0.  ) 

YBLO  = 

YR ( 2 )  ♦  DB 

/ 

2 

0 

.  AND  . 

YR  ( 2  ) 

.GT 

0.  ) 

YBLO  = 

YR ( 2 )  -  DB 

/ 

2 

0 

.  AND  . 

YR  ( 2  ) 

.LT  ‘ 

0.  ) 

YBLO  = 

YR ( 2 )  +  DB 

/ 

2 

C 

C 


C 

c 

c 


160  START  =  .FALSE . 

WR I T  E ( 6 , 7  50 )  XR 

170  FORMAT  (2*,  'REVERSING  POINT  AT  XR= ' .  E16.5/) 

WRITE (6 . 400 )X .MODE ( 1 ) .MODE (2),YP(1).YP(2).YBLI,YBL0 
K  (  J 1 1  )  =  MODEDB  -  1 
YBLOS  =  YBLO 
RETURN 


DETERMINE  THE  POINTS  OF  LEAVING  THE  DEAD-BAND 


180 

IF 

(  .NOT. 

START) 

GO  TO 

190 

I  F 

( DAE  S (YBLI)  .GT 

.  DE/2 

.  )  GO 

GO 

TO  220 

190 

I  F 

(J1  .LT. 

0)  GO 

TO  140 

IF 

( MODE ( 1 ) 

.  EO.  1 

.  AND  . 

YBL  I 

IF 

(MODE (  1  ) 

.  EO.  - 

1  .AND.  YB 

GO 

TO  220 

TO  200 


.GT.  ( YR ( 2 )  +  DB))  GO  TO  200 
I  .LT.  ( YR ( 2 )  -  DB))  GO  TO  200 


200  MODEDE  =  0 

CALL  DDB AND ( Y .  YP.  YS,  YSP.  X.  H.  YR(2),  XDB.  DB .  N.  M.  START) 
C  WR I TE ( 6 , 2  3 ) XDB 

210  FORMAT  ( 2X ,  'DEAD-BAND  LEAVING  POINT=',  E16.4/) 


a  ) 
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197 

198 

199 
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201 
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203 
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205 
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2  10 

2  1  1 

2  1  2 

213 

214 

2  15 

2  16 

217 

2  18 

2  19 

220 

22  1 

222 
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227 

228 

229 

230 

23  1 

232 

233 

234 

235 

236 

237 

238 

239 

240 


n  n  o 
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c 

c 

c 


c 


c 


CHECK  FOR  REFLECTION  POINT  JUST  AFTER  LEAVING  THE  DEAD-BAND 
STATUS  =2 


MODE (  1  )  =  MODE ( 2 ) 

IF  (YSP(M)  .GT.  0)  MODE ( 2 )  =  1 
IF  (YSP(M)  .LT.  O)  MODE ( 2 )  =  -1 
J1  =  MODE (  1 )  *  MODE ( 2  ) 

IF  (J1  .LT.  0)  GO  TO  140 
YP ( 1 )  =  YP ( 2  ) 


YP ( 2 )  =  YSP(M) 

Y( 1)  =  Y ( 2 ) 

Y ( 2 )  =  YS(M) 

IF  ( V  s ( M )  .LT.  0  .AND. 
IF  (YS(M)  .GT.  0  .AND. 
IF  (YS(M)  .LT.  0  .AND. 
IF  ( Y  S ( M )  .GT.  0  .AND. 
WR I T E ( 6 , 400 ) X . MODE (  1  )  , 
K( J1 1 )  =  MODEDE  -  1 
YBLOS  =  YBLO 


MODE ( 2 )  .GT.  0)  YBLO  =  YS(M)  +  DB  / 

MODE ( 2 )  .LT.  0)  YBLO  *  YS(M)  -  DB  / 

MODE ( 2 )  .LT.  0)  YBLO  =  YS(M)  +  DB  / 

MODE ( 2 )  .GT.  0)  YBLO  =  YS(M)  -  DB  / 

MODE ( 2 ) . YP( 1) ,YP(2) , YBL I .YBLO 


2  . 
2  . 
2  . 
2  . 


RETURN 

C 

220  Y ( 1 )  =  Y ( 2 ) 

Y ( 2 )  =  YBL  I 
YP (  1 )  =  YP ( 2  ) 

\P<2)  =  YPBLI 

C  WR I T  E ( 6 . 400 1 X , MODE ( 1 ) , MODE (2),YP(1).YP(2),YBLI  YBLO 

230  DO  240  I  =  1  ,  N 

YSP( I  )  =  Y I P ( I  ) 

240  YS< I )  =  YI (I) 


STATUS  =  1 
K  (  J 1 1  )  =  MODEDE  -  1 
YBLOS  =  YELO 
RETURN 


CHANGE  IN  THE  INPUT  WITIN  THE  DEAD-BAND 


250  YBLO  =  YBLOS 
MODEDE  =  1 
C  WR I TE ( 6 , 200 ) 

260  FORMAT  (/'CHANGE  IN  THE  INPUT  WITHIN  THE  DEAD-BAND'/) 

C  WRITE (6 ,400)X .MODE ( 1 ) .MODE (2 ) , YP(  1  )  ,  YP ( 2  )  . YBL I . YBLO 

K( J1  1  >  =  MODEDE  -  1 
RETURN 
END 
C 

SUBROUTINE  DDB AND ( Y .  YP.  YS,  FS,  X,  H.  YR.  XDB .  DB .  N.  M.  START) 
IMPLICIT  REAL  *8 ( A  -  H.O  -  Z) 

REAL  *8  Y ( 2  )  .  YP(2).  YS(N),  FS(N) 

INTEGER  RUNGE1 

COMMON  /BL0CK3/  YR 1 ,  M0DE(2).  MODEDB 
LOGICAL  START 
C 

YR 1  =  YR 

XI  =  X  -  2 .  *  H 

X2  =  X  -  H 

IF  (  .NOT.  START)  GO  TO  10 

IF  ( Y ( 2 )  .LE.  0.0)  CO  =  YR  -  DB  /  2. 


. 


24  1 

242 

243 

244 

245 

246 

247 

248 

249 

250 

25  1 

252 

253 

254 

255 

256 

257 

258 

259 

260 

26  1 

262 

263 

264 

265 

266 

257 

268 

269 

270 

27  1 

272 

273 

274 

275 

276 

277 

278 

279 

280 

28  1 

282 

283 

284 

285 

286 

287 

288 

289 

290 

29  1 

292 

293 

294 

295 

296 

297 

298 

299 

300 


non  o  o  n  o 
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c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 


10 


IF  (  Y  (  2  )  .GE. 
START  =  .FALSE 
GO  TO  20 
IF  (YR  .LT.  0. 
IF  (YR  .LT.  0. 
IF  (YR  .GT.  0. 
IF  (YR  .GT.  0. 


.0)  CO  *  YR  ♦ 


.AND.  MODE ( 2 ) 
•AND.  MODE (2) 
.AND.  MODE ( 2 ) 
.AND.  MODE ( 2 ) 


DB  / 

2  . 

.GT. 

0) 

.LT. 

0) 

•  GT. 

0) 

.LT. 

0) 

CO 

CO 

CO 

CO 


YR 

+ 

DB 

YR 

- 

DB 

YR 

+ 

DB 

YR 

- 

DB 

PERFORM  INVERSE  INTERPOLATION  ;T0  GET  THE  TIME  AT  WHICH 
THE  INPUT  IS  LEAVING  THE  DEAD  BAND. USING  NEWTONS  FORMULA 
(SECOND-ORDER) 


20  XDB  =  XI  +  . 1  *  H 
J  =  1 

30  A  =  2  '  (  Y  (  1 )  -  Y  ( 2 ) )  /  ( X2  -  XI)  **  3  +  (YP(2)  +  YP (  1 ) )  /  (X2  - 
1X1)  **  2 

B  =  (3  *(X2  +  X 1 ) *  (  Y  ( 2 )  -  Y (  1  )  )  )  /  (X2  -  XI)  **  3  -  ((2.*X2  +  XI)* 
1  Y  P (  1  )  +  ( 2  .  *  X  1  +  X2  )  *  YP ( 2 ) )  /  ( X2  -  XI)  •*  2 

C  =  (6.*X2*X1*(Y(1)  -  Y ( 2 )  )  )  /  ( X2  -  XI)  •*  3  +  ((X2**2  2.*X2* 

1 X  1  )  *  YP (  1)  +  ( X 1  *  *  2  «■  2  .  *X2*X 1 ) * YP( 2 ) )  /  (X2  -  XI)  **  2 
D  =  ( ( X 2  -  3  .  *X  1  )  »X2**2* Y( 1 )  +  ( 3  .  * X2  -  X  1  )  *X  1  * *2 • Y ( 2 ) )  /  ( X 2  - 

1X1)  **  3  -  (X2*X1*(X2*YP(  1)  X 1  *  YP  ( 2 )  )  )  /  (X2  -  XI)  **  2  -  CO 

SOLVE  A/*»3+BX**2+CX+D  USING  NEWTON'S  METHOD 
THE  ANALYTICAL  SOLUTION  OF  THE  OUBIC  EQUATION  CAN  BE 
DEDUCED  BUT  THE  ACCURACY  OF  ITS  NUMERICAL  SOLUTION  IS 
VERV  POOR. 

K  =  1 

40  FSS  =  A  *  XDB  *  *  3  *  B  •  XDB  *  *  2  *  C  *  XDE  +  D 

DFSS  =  3  DO  *  A  *  XDB  *  *  2  ♦  2 . DO  •  B  *  XDB  +  C 

IF  (DAES(FSS)  -  ID-5)  60.  60.  50 

50  XDB  =  XDB  -  FSS  /  DFSS 
K  *  K  +  1 

IF  (K  .GT  30)  GO  TO  60 
GO  TO  40 

60  IF  (XDE  LT.  X  .OR.  XDB  .GT.  X2)  GO  TO  80 
WRITE  (6,70) 

70  FORMAT  ('THERE  IS  NO  ROOT  IN  DDE  AND ' / ) 

STOP 


DETERMINE  BACK  LASH  LEAVING  POINT 
INTEGERAT  USING  RUNGE-KUTTA  METOHED . 


80  HI  =  XDB  -  X2 
XRUNGE  *  X2 

90  K 1  =  RUNGE 1 (N. YS . FS . XRUNGE .HI ) 
IF  ( K 1  .  NE  1 )  GO  TO  100 
CALL  DE1KYS.  FS.  XRUNGE) 

GO  TO  90 

100  CALL  DE  1  1  ( YS.  FS.  XDB) 


CONTINUE  INTEGRATION  FROM  XDB  TO  X 


MODEDB  *  0 
HI  =  X  -  XDB 


301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

31  1 

312 

313 

3  14 

3  15 

3  16 

317 

3  18 

3  1  9 

320 

32  1 

322 

323 

324 

325 

326 

327 

328 

329 

330 

31  1 

332 

333 

334 

330 

336 

337 

338 

339 

340 

34  1 

342 

343 

344 

345 

346 

347 

348 

349 

350 

35  1 

352 

353 

354 

355 

356 

357 

358 

359 

360 


o  o  o  o  o  o 
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110  K 1  *  RUNGE 1 (N, YS , FS . XDB . HI ) 

IF  ( K 1  .NE.  1)  GO  TO  120 
CALL  DEI  1 ( YS ,  FS.  XDB ) 

GO  TO  1  10 

120  CALL  DE  1  1 ( YS ,  FS.  X) 

RETURN 

END 

SUBROUTINE  I NTR ( Y ,  YP.  YS.  FS.  X.  H.  YR .  XDB.  XR.  N  M) 
IMPLICIT  REAL  *8 ( A  -  H.O  -  Z) 

REAL  *8  Y ( 2  )  .  YP ( 2  )  .  YS(N),  FS(N) 

INTEGER  RUNGE1 

COMMON  /BL0CK3/  YR 1 ,  M0DE(2).  MODEDB 


FIND  THE  REVERSING  TIME  BY  IMPLEMENTING  HERMITE  INVERSE 
INTERPOLATION  FORMULA  ( OUB I C- INTERPOL AT  I  ON ) 


C 


c 


c 


C 


c 


c 

c 

c 


c 

c 


J  =  1 
YR  1  =  YR 


X  1 

=  X  - 

2  .  * 

H 

X  2 

-  X  - 

H 

XR 

=  (XI 

+  X2) 

/  2. 

10 

A  = 

( -6 

*  (  Y  (  2  ) 

-  Y(  1  )  ) 

♦  3 

1  *  * 

3 

B  = 

(6  » 

(  X  2  -*■ 

X 1 ) * ( Y( 2  ) 

1  ♦ 

( 2  -X2  ♦  XI 

)  *  Y  P  (  1  )  )  ) 

/ 

C  = 

(  -6 

*X2*X  1 

* ( Y( 2)  - 

Y  (  1 

1  + 

(  X2*  * 

2*2. 

•X2*X  1  )  » YP(  1 

K  = 

1 

20 

FSS 

=  A 

’  XR  * 

*  2  *  B  • 

XR 

DFSS  =  2 

.  *  A 

»  XR  «■  B 

IF 

( DABS ( FSS ) 

-  ID-5) 

40. 

30 

XR 

-  XR 

-  FSS 

/  DFSS 

K  =  K  *  1 

IF  (K  GT.  30 1  GO  TO  40 


40  IF  (XR  .LT.  X  OR.  XR  . GT .  X2)  GO  TO  60 
WRITE  (6.50) 

50  FORMAT  ('THERE  IS  NO  ROOT  IN  I NTR ' / ) 
STOP 


INTEGRATE  THE  SYSTEM'S  DYNAMICS  FROM  X2  TO  XR 

60  HI  *  XR  -  X2 
XRUNGE  =  X2 

70  K  1  =  RUNGEMN.YS.FS,  XRUNGE.  HI) 

IF  (K 1  .NE .  1 )  GO  TO  80 

CALL  DE 1 1 ( YS ,  FS.  XRUNGE) 

GO  TO  70 

80  CALL  DE 1 1 ( YS .  FS,  XR) 

X2  =  XR 
YR  *  YS(M) 


CONTNUE  INTEGRATION  FRO  XR  TO  X 


■ 


36  1 

362 

363 

364 

365 

366 

367 

368 

369 

370 

37  1 

372 
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c 

STATUS  =  1  . 

MODEDB  =  1 
HI  «  X  -  XR 

90  K 1  *  RUNGE 1(N,YS.FS.XR,H1) 
IF  (K1  .NE.  1)  GO  TO  100 
CALL  DE1KYS,  FS.  XR) 

GO  TO  90 

100  CALL  DEI 1( YS.  FS.  X) 

C 

RETURN 

END 


■ 
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